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PRECEDING PTOB BEX** m r 


ANALYSIS OF ESTIMATION ALGORITHMS 
FOR CDTI AND CAS APPLICATIONS 

Tsuyoshi Goka 

Analytical Mechanics Associates, Inc. 

Mountain View, California 94043 

SUMMARY 

The objectives of this project were to analyze and/or to develop 
estimation algorithms for Cockpit Display of Traffic Information (CDTI) 
and Collision Avoidance System (CAS) applications. The algorithms are 
based on actual or projected operational and performance characteristics 
of an Enhanced TCAS II traffic sensor developed by Bendix and the Fed- 
eral Aviation Administration. 

Three algorithms are examined and discussed. These are horizontal 
x and y, range and altitude estimation althorithms. Raw estimation er- 
rors are quancified using Monte Carlo simulations developed for each 
application; the raw errors are then used to infer impacts on the CDTI 
and CAS applications. Applications of smoothing algorithms to CDTI pro- 
blems tire also discussed briefly. 

Conclusions are summarized based on the analysis of simulation 
results. 
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INTRODUCTION 

NASA Langley Research Center is pursuing a research effort con- 
cerning the Cockpit Display of Traffic Information (CDTI) concept [1). 
The CDTI is a device which presents information to the pilot and crew 
depicting the status of surrounding traffic including position and 
velocity states. The traffic information is provided by a "traffic 
sensor." Most promising candidate sensors are FAA developed Traffic 
Alert and Collision Avoidance Systems (TCAS) . 

TCAS is strictly an airborne system which provides the aircraft 
separation protection information independent of the ground ATC sys- 
tem. The FAA plans call for developing two types of TCAS— TCAS I and 
TCAS II. Within each category, a certain latitude in capability is 
allowed to satisfy a wide spectrum of user requirements. The enhanced 
TCAS II which is capable of obtaining relative bearing measurements 
may be able to support CDTI applications. There are two designs in 
this enhanced TCAS II category. One design developed by MIT/Dalmo 
Victor is based on the so-called active Beacon Collision Avoidance 
System (BCAS) . The other developed by Bendix is based on the so- 
called full BCAS concept. The former unit is being, tested in actual 
commercial flight operation environments; and the other is undergoing 
an extensive flight test with the prototype system. 

The current effort is a part of parallel efforts consisting of: 

(a) Development of a realistic enhanced TCAS II simulation 
model, and; 

(b) Analysis of the TCAS estimation algorithms for the CDTI 
and CAS applications. 

The companion report, "Enhanced TCAS IX/CDTI Traffic Sensor Digital 
Simulation Model and Program Description [2], contains a detailed dis- 
cussion of the Bendix designed system. A shorter version of TCAS II 
description is given In Appendix A. 
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The objectives of the current effort are to analyze and/or develop 
estimation algor ithns for CDTI and possibly for Collision Avoidance 
System (CAS) logic applications based on actual or projected operational 
and performance characteristics. For purposes of estimation three coordinate 
axes are considered separately. These are x and y, range and altitude. 
Within these axes, variations in signal configurations, signal sources 
or algorithm implementation items are examined. 

The following procedure is used to compare and analyze the esti- 
mation algorithms. The raw estimation errors are quantified using Monte 
Carlo simulation method. The raw error data are then used to infer im- 
pacts on the CDTI and CAS applications. For example, altitude and alti- 
tude rate estimation errors per se do not mean much; however, if these 
are factored into the projected altitude error, then the latter would have 
a great significance in terms of safe altitude separation threshold of CAS 
logic . 

In Chapter II, basic performance of horizontal x-y filter algorithms 
are discussed. The basic question in this chapter is how and what kinds 
of other signals (in addition to the relative range and bearing measure- 
ments) are best utilized to provide better estimates in horizontal x-y 
axes. This is motivated by the fact that the dynamic lag due to maneuvers 
by Own or target aircraft induces large and sustained errors in position 
and velocity estimates. These errors can be compensated by utilizing 
maneuver parameters (such as heading angles) in estimation algorithm. This 
assumes that the target data are made available via the Mode S data link 
capability. Also, the questions of filter gains determination are ad- 
dressed. The filter gains depend on many operational factors; thus, it 
is not o trivial matter. These questions are probed by means of error 
statistics generated by Monte Carlo simulation program. 

In Chapter III, range and range rate estimates are obtained in 
several ways. Raw error statistics for each are obtained by Monte Carlo 
method using "realistic" encounter scenarios. These are, in turn, 
compared and analyzed in terms of accuracy. 
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In Chapter IV, an altitude tracker algorithm is developed and 
presented. The altitude axis poses a special estimation challenge in that 
the target altitude measurements are quantized to the nearest 100 ft. 

This causes a certain observability problem. The algorithm is based on 
the level switching time detection concept. The performance analysis are 
carried out by comparing the estimation errors with those of a non-quantized 
alpha-beta tracker. The latter represents the best possible without aiding 
the estimation algorithm with external signals. 

In Chapter V, the raw error statistics (obtained in the previous 
three chapters) are analyzed from the user's view-point, i.e., from the 
CDTI and CAS logic application aspects. These would provide relative and 
absolute merits of particular estimation algorithms with respect to the 
design requirements. Also, a short discussion of smoothing (rather than 
estimation) algorithms are given with respect to CDTI applications. 

Appendices A through D provide peripheral but important information 
which are directly related to this effort. 
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II 

HORIZONTAL X-Y FILTERS 
Filter Configuration 

In the context of an airborne CDTI sensor based on the proposed 
enhanced TCAS II, the relative position of an intruder aircraft is obtained 
in range and bearing axes. The verticle axis is provided by an encoding 
altimeter (either above MSL or AGL) . The range and bearing measurements 
are taken with respect to a cylindrical coordinate system attached to Own 
fuselage; therefore, the coordinate system is itself subject to transla- 
tional and rotational motions as Own aircraft undergoes maneuvers. 

The sensor is designed to account for Own's orientation effect by 
means of software compensation utilizing onboard INS generated attitude 
angles. (A brief description of the sensor surveillance and operational 
characteristics is given in Appendix A.) For filter analysis purposes, 
point-mass kinematics are assumed for both Own and Target aircraft. That 
is, that the effect of Gum's orientation angles is assumed to be negligible. 
This may be justified by the fact that (a) an INS provides accurate orienta- 
tion, and (b) low frequency bias errors do not affect velocity estimates. 

If the altitudes are ignored, the relative horizontal measurements (range 
and bearing) are given by 

r m - < Ax2 + A y 2]l5+ ^ r . (1) 

b m - tan -1 (Ay/Ax) + $ b , 

where £ r and £ b are measurement errors. 

For the purpose of designing estimation algorithms, a model is needed 
to describe the relative kinematics. Now, the magnitude of purely longi- 
tudinal acceleration, (i.e., along the velocity vector) is small - usually 
less than 2-3 kt/sec for commercial operation. Thus, the longitudinal 
acceleration effect (approximately 2.5 ft in position change at 3 kt/sec) 
is masked by somewhat large measurement errors (ranging error of 100 feet). 
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However, the acceleration due to a turn is not so trivial. 


A 1/2 g acceleration corresponding to the bank angle of 25 deg is not 
unconmon for commercial operation. A suitable kinematic model is obtained 
by assuming that aircraft follow a series of straight line or circular arc 
segments. If position and velocity vectors, p and v, are defined with 
respect to an earth fixed rectangular coordinate system, then each aircraft 
is described by an equation of the form 

i 


*[•*}[• ‘‘iw 


i - (0) Own or Target (T) 


( 2 ) 


-*> 


t ”o] - anlI ’P >] • 


The turn rate, is a piece-wise constant time function. By sub- 
tracting, the relative kinematic equation for two aircraft can be expressed 


as 


d 

dt 


Ap 


0 I 


Ap 

+ 

0 

- Av 


1 

O 

o 

1 


Av_ 


Aa 


(3) 


where Aa is the relative acceleration given by 


(4) 


Aa - n T v T - n 0 v 0 - + (fl T - n 0 )v o . 

Obviously when both aircraft are non-accelerating (straight line flight), 
then Aa ■ 0, or Av - constant. 


Equations (1) and (3) form the basis for designing horizontal x-y 
filters to estimate position and velocity. In the following sections 
several filter algorithms are derived and discussed. The configuration 
differences are based on different signals available to aid the relative 
position measurement. These include: 
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(1) no additional measurements are available; 

(2) Own body rates are available; and 

(3) Own and Target body rates are available. 

In cases (2) and (3) it is assumed that Own or Target acceleration 
can be derived to compensate for the unknown acceleration components. In 
case (3), availability of suitable avionics and an air-to-air communication 
link is assumed. Because a Mode S transponder is an Integral element of a 

TCAS t it is only necessary to establish a digital data link between an onboard 
sensor such as an INS or a navigation computer to the so-callod Standard 

Message Interface port of the Mode S transponder. It is noted that the 
Own signal compensation is relatively easy in the sense that a digital 
data interface between an INS and TCAS is in place. Cross-feed signals 
include the three Euler angles <|>, 6 and i|> and the ground speed, . 

Own acceleration can be generated easily from these available signals. 

Non-Aided Filter Configuration In the absence of any acceleration 
indication, there are essentially two approaches for designing estimatio'i 
algorithms. One is simply to ignore the acceleration input. The other 
is to estimate the unknown acceleration. The second approach is generally 
very difficult. Reference [3,4 ] discusses several methods applied to a 
simpler problem of tracking maneuvering aircraft using a ground based 
Mode S sensor. Compared to the previous ground based study, there are major 
differences in TCAS surveillance functions : 

(1) ground-based versus aircraft based ; 

(2) sampling period of 4-5 sec vs 1-8 sec; and, 

(3) bearing error of 0.04 deg vs 1-2 deg. 

The most crucial airborne disadvantage is the bearing error which is 
25-50 times larger in magnitude. The linear equivalance is 65 ft compared 
to 1600-3200 ft at a range of 15 rani. Therefore, attempting to estimate 
unknown acceleration in the given noise environment is not realistic. Thus, 
the first approach is now studied further. 
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If the acceleration is assumed to be zero, then Eq (3) reduced to 


_d 

dt 


AP 

Av 


:][::] 


( 5 ) 


By discretizing over one sample time interval, At, and by writing 
the x component equation (y equation is identical), then it follows 


Ax 

Ax 


n+1 


i I Tax 
1 J U* J n 


( 6 ) 


The pesudo measurement equations are 

*x » 

Using Eqs. (6) and (7), the standard filter algorithm is given by 


A ® 

Ax_ - 

m 

r 

cosb m = 

Ax^ 

+ 

n 

n 

n 

n 


< - 

m 

r 

n 

sinb m - 
n 

Ay 

+ 


(7) 


position prediction: 

Ax + 

' : 

■ Ax + At * Ax ; 

n n 


measurement error: 

Ai 

* Ax n+1 ~ Ax+ ; and 

(8) 

estimate update: 

Ax * 

Ax + + 8xi • 



* 

Ax • 

Ax + 8 x 2 *^ x 



where g x ^ and g x 2 are filter gains. Note that Eq. (8) would represent an 
optimal filter if Own and Target aircraft are not maneuvering. Choosing 
proper filter gains is more of an art than a science. There are many 
factors involved such as noise variance reduction, dynamic error minimiza- 
tion, computational ease, and so on. This is discussed in a later section. 


Own Signal-Aided Filter Configuration When the Own acceleration 
signal is available, at least half of the acceleration term in the dynamic 
equation can be compensated. Thus, after proper discretisation, the model 
equation equivalent to (6) is given by 
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( 9 ) 
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At A 

1 A 


Ax 

Ax 


At 2/2 ] * 

At J °» n 


where x 


0,n 


is the estimated acceleration from Own onboard signals. 


Taking advantage of an inertial grade INS, differential ground velocity. 


£v x can be used. 


6v_ 


« 

x. 


- x„ 


* At 


'x,n ' v 0,n 0,n-l 

The estimation algorithm is given by 


0,n 


( 10 ) 


state prediction: 

> 

*+ 

i 

Ax_ + 
n 

At . Ax_ 
n 

- (At/2) 6v , 

A j H 



Ax + * 

> 

> 

1 

6v ; 

x,n * 



measurement error: 

Ax ■ 

A m 

Ax n+1 

. + 

- Ax ; 

and 

(id 

estimate update: 

Ax . , 
n+1 

■ Ax m , . 
n+1 

+ 8x1 

• Ax 9 



4 Vu 

- Ax + 

+ 8x2 * 

Ax , 



where g x ^ and g x 2 are again the filter gains. As expected, this filter 
configuration would be ideally suited if the Target aircraft does not 
accelerate. Also, an intuitive expectation that this filter would be 
"twice as good" as the previous configuration because it compensates 
for at least one half of the problem causes is not correct. Even though 
this is true most of the time, it is apparent from Eqs (3) and (4) that 
if the relative acceleration Aa - 0 even though Sq and are large, then 
the Non-Aided Configuration would outperform the Own-Aided Configuration. 


Own and Target Signal-Aided Filter Configuration When Own and Target 
signals are both available, then most of the unknowns in the system dis- 
appear. For the purpose of algorithm design, it is assumed that the 
Target cross-links its ground speed and heading ^ as part of its 
Mode S surveillance reply message. (If the Target is equipped with an 
enhanced TCAS II, then these signals are available per specification, i.e., 
no hardware modification is necessary.) Under these assumptions, Target 
differential ground velocity can be computed as 
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COS(* T „) 

T,n 
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6v Tx,n " v GT,n 


- v 


GI,n-l co,< *T,n-l ) * 11 ' *Tn 

( 12 ) 


6v Ty,n " v GT,n #1 ®<*t,ii* " V GT,n-l “^^T.n-D “At * y^ 


It is noted that the sampling time for Own signals and that of Target 
signals may not be the same, i.e., they are not synchronous. Target data 
are available only when the interrogation/repiy cycle is complete, which 
may take from 1 to 8 sec; Own signals are available at the TCAS basic 
cycle time of 1 sec. For example. Target data may be available every A 
sec compared to Own data available every 1 sec. This Implies that special 
care is needed in processing Own inertial data. 


If the computed acceleration terms are incorporated, the discrete 
system equation becomes 


Ax 


1 At 


Ax 


At/2 


At/2 
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The corresponding estimation algorithm is given by 


state prediction: 

Ax + - 

Ax + 
n 

At • 
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ix f n 
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measurement error: 

- A Vi • 
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9 
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(14) 

estimate update: 
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where g . and g x 2 are filter gains. 



As far as the estimation algorithms are concerned, there are very 
small computational differences. The essential difference is in computing 
the predicted state; the rest is identical. The real computational loads 
for the aided configurations lie in pre-processing and keeping track (time- 
and book-keeping) of Own and Target Inertial signals. These seem rather 
trivial compared to other more complicated processing performed in other 

TCAS modules. 


Filter Gain Generation 

Once the configuration is chosen, the filter gains need to be specified 
in order to implement the algorithm. In this section, three methods of sel 
ecting the gain values are discussed. These are (a) fixed 08 tracker gain, 

(b) Kalman filter gain, and (c) table-look-up gain. These methods are ana- 
lysed for each of the three configurations. 

Method (a) would be simplest from the computational point of view. The 
fixed gain configuration would have a major disadvantage: the fixed filter 

gains imply the noise reduction ratio remains the same. Thus, the "roughness" 
of the estimate would be proportional to the "roughness" (or noise magnitude) 
of the input signal. As is well-known, the noise variance of the pseudo x or 
y measurement is affected by the so-called range effect due to the basic 
radar coordinate system. Thus, the input variance would be proportional to 
range; hence, the estimate error variances would be proportional to range. 

Method (b) would be most suitable in terms of best performance in the 
sense that gains are automatically adjusted according to input noise variance. 
This feature would "desensitize" the range effect problem encountered by the 
first method. The price for the added performance gain is additional compu- 
tational load. 

Method (c) tries to strike a mid-point in performance and in computa- 
tional load. Very briefly, the filter time constant (approximately the 
reciprocal of the filter bandwidth) parameter is stored in a two-parameter 
look-up-table of the sampling period and measurement noise standard devia- 
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tion level, 
back gains. 


The tine constant is used to compute position and velocity feed- 
For convenience, important equations are repeated in Table 1. 


A few points need to be noted. The estimated acceleration, ® n in 
Eq (15) is 0 for non-aided configuration; Own estimated acceleration, 

A A A 

a^ for Own-aided configuration; or a Tn - a^ for Own-and-Target-aided 
configuration. The model dynamics and state prediction Eqs (15) and (17) 
are identical; therefore, the modeling error is dominated by the accelera- 
tion error. 


As pointed out previously, the pseudo measurement errors ^ 
Eq (16) would have the covariance matrix 


( 21 ) 

(and E ) in 
yn 






2 

o o 

x xy 

o a 2 

L *y y J 

" o 2 2 2 2 

cos b o + r sin b a u 
r d 

2 2 2 
cosb sinb (o - r a, ) 
r d 


( 22 ) 

2 2 2 

cosb sinb o_ -r o. 

r D 

2 2 2 2 
sin b + r cos b o 



Here o r and are range and bearing error standard deviations. Note that 
the off-diagonal term is generally non-zero. 


Method (a) - Fixed Gains This is the method used to obtain the a and 
B tracker gains in the current TCAS design. Therefore, the method is directed 
toward the non-aided configuration; however, the design procedure can be 
applied to other configurations. The following discussion follows Ref. [5] 
very closely. The basic idea is to compute the error due to measurement 
error and the error due to acceleration error. If these errors are com- 
bined statistically in a correct way, then we can optimize the gain values 
to minimize the total error. 


(i) Error due to measurement error only 
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Estimation Model 




Table 1. Summary of Estimation Equations 


Dynamic Eqn. 


Ax 1 Lt Lx Lt 2 /2 

m # + 

Lx. n+1 .0 1 n L 


Measurement Eqn. 


m . m 
r cosb^ 
n 


Ay - r sinb, 
n n 


- Ax + £ 

n x,n 


S Ay + C 

n y.n 


State Prediction 


Measurement 

Error 


. m . + 
Ax - Ax 
n 


State Update 


Estimation Error 


_ Ax"| fAxl + r a 1 
_A*x n+1 L A *2 : J LsMtJ 



AxJ n+1 


1-a At(l-a) | Ax 

_-6/At 1-6 J _Ax_ n 


a 

+ t 

6/At x »' 


** Ax 


■ 0 (Non-aided) ; “ -a 0x n (O wn “ a lded) ; and 

- *T*,„ - 4 0x,n <°“" “ d 

= Ax - Ax ; Ax = Ax - Ax (state estimation error) 


a - a 


(acceleration model error) 


If the acceleration error, a R in Eq (21) is set to zero, and we solve 

for the steady-state error covariance using the standard linear covariance 

2 

propagation foraula (assuming o x is stationary), then the following re- 
lationships are obtained: 



o 


:2 

X 


2a 2 + B(2-3o) 1 2 

a(4-B-2a) J a x * 



(23) 


These two formulas express the estimation errors in terms of input error 
variance, sampling period and the filter gains • 

(ii) Error due to acceleration model error. 

Now, set the input noise to zero. Then, Eq (21) expresses the esti- 
mation errors due only to the acceleration error. Since the nature of this 
error is low frequency, we can compute "worst" type error due to "worst" 
possible acceleration error. By solving for the steady state error, given a^ 
as the maximum model error, the following relationships are obtained ; 


x 


ss 


7 

X 


ss 


(At ) 2 - 

e * 

At 2a-B - 
2 B ®M ’ 


(24) 


These two formulas express the steady state errors in terms of maximum 
/ acceleration uncertainty, sampling period and the filter gains. 


(iii) Total Errors 

A natural way to combine the statistical errors of Eq (23) and the 
deterministic errors of Eq (24) is to take the root mean square. The 
individual errors can be interpreted as Eq (24) being the mean error 
and Eq (23) being the variance about the mean. Thus 
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a(4-B-2a) 


4 (l-ar 


rms | A fc 2 a(4-6-2a) x 


Formulas (25) express the to:al estimation errors in terms of measurement 
noise m. >nitude (o^) , maximum acceleration uncertainty (a^) , the sampling 
period (At) and the filter gains (o and 6) . Therefore, given suitable 
numbers for o x , a^ and At, the best a and 8 can be obtained which minimizes 
the rms errors. In Ref. I 5 ], this was done by using o x « 825 ft, and 
^ - 0.5g “ 16.1 fps 2 for At - 1 ',...,8 sec. Table 2 shows the "optimum” 
gains. 

Table 2. Optimum Gains for Non -Aided Configuration 


At (sec) 




8/At 

-1 


O'- (ft) 406 



*„ <“> 1 183 



X (ft) 445 

rms 1 



2 

3 

0.37 

0.465 

0.0875 

0.1 

506 

570 

95.5 

102.8 

1 

232 

258 

70.2 

50.7 

557 

626 

118.5 

114.6 





5 

6 

0.58 

0.62 

0.113 

0.114 

650 

676 

114,8 

117.4 

299 

322 

42.4 

39.1 

! 

716 

749 

122.4 

123.8 
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The value o£ o - 825 ft was chosen, because this number has been "traditionally" 
important in the airborne based collision avoidance concept 16 ]. The maximum 
acceleration a„ * 0.5g was selected based on the maximum appearing in the so- 
called MITRE's FAA basic model 1 7 ). Some comments are needed. 

A ° x of 825 ft implies that the target range is 7.8 nmi if bearing error 
is 1 deg and 3.9 nmi if 2 deg. Beyond these ranges, the actual linear error 
would be larger (e.g., at 10 nmi range and 2 deg bearing error, o x « 2100 ft). 
Within these ranges the linear error would be smaller (e.g., at 2 nmi range 
and 1 deg bearing error, * 210 ft). Therefore, the range effect needs to 
be accounted for. 

In the steady state error derivation, it was assumed that the maximum 
acceleration was maintained indefinitely. This assumption represents the 
worst possible case. From Table 2 it is clear that the position and velocity 
errors are dominated by the measurement noise magnitude but not by the 
dynamic error. This seems to indicate that the gain values could be made 
a little lower at this measurement error level. However, it should be noted 
that the above gains were derived with Own being stationary. Thus, if Own 
aircraft is also conducting a turning maneuver, the acceleration uncertainty 
could be larger. 

Similar analysis may be pursued for other configurations. But it 
would be more effective to proceed to the other two gain selection methods. 

Method (b) - Kalman Filter Gains One of the major disadvantages of 
a fixed gain filter in a radar environment is that the range effect is not 
automatically accounted for. That is, "optimized" gains at the error level of 
a x - 825 ft (range of 3.9 nmi at Ojj - 2 deg) is optimal at that point and 
suboptimum everywhere else. This is the main motivation for utilizing gains 
based on the Kalman filter theory. 

There are two difficulties with this approach which need to be discussed. 
One is the treatment of acceleration uncertainty. The other is the coupling 
problem of x and y axes due to non-zero covariance (o x y t 0) in the pseudo- 
measurement errors. The latter problem can be significant to the extent that 
the state covariance should be a 4 x 4 rather than two 2x2 matrix. 
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For the sake of saving computations, this statistical dependence Is 
ignored, and the filter gains are solved only fox decoupled axis. 


The problem of unknown acceleration is not trivial. In the usual 
application, the process noise magnitude is varied to "tune" the Kalman 
filter. Therefore, the unknown acceleration is assumed to be a white noise 
process, and the magnitude is varied so that the gains yield a satisfactory 
filter performance. For the purpose of obtaining the Kalmar filter equation, 
the acceleration input in Eq (15) is assumed to be a zero-mean white noise 
process with standard deviation of o a . (Equation (17) needs to be modified 
accordingly.) Furthermore, the measurement error is approximated by 

o 2 = r 2 sin 2 b a 2 ; a 2 = r 2 cos 2 b o 2 . (26) 

x by b 

With these simplifying assumptions, the usual covariance equations for 
line? . systems can be developed. For the two state filter, the equations 
are particularly simple [ 8 ]. 


Covariance Prediction (Propagation): 


+ 

p l 

= Pl + 2At p 2 + At 2 

P 3 + 

(At 3 /3 ) o a 2 » 


+ 

p 2 

P 2 + At 

P 3 + 

(At 2 /2)o a 2 , 

(27) 

4 



2 


p 3 

* 

P 3 + 

At o„ 

a 



Gain Computation: 

4 - , + 2.-1 

“n ' Pj <Pj + % ) 

2 + 2-1 
B n /At - p 2 ( Pl + a x ) 


(28) 
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Covariance Update: 

2 

p, » a o 

F 1 n x 

2 

p 2 - CS n /At)a x 

p 3 * P 3 ” <e n /^t ) 2 3 <p* + a x 2 ) 


(29) 


The terms are defined as 

, *2 

Pl = E(Ax 2 ) , p 2 * E(Ax Ax) , and p 3 = E'Ax ) , 

and (.)"* indicates the predicted value according to the usual linear 
system propagation formula. Time reference is suppressed in the above; 
however, if the equations are coded in the exact sequence appearing above 
then the recursive nature of algorithm is maintained. Furthermore, the 
Pi +, s can be stored in the p^'s locations. 


The conclusion is clear. These computations are simple enough that 
implementation in a micro-processor should be straight forward. In an 
actual Implementation, some other considerations need to be made. 


2 2 

(1) Equation (26) cannot be used to compute o x and o y because 

the true values of r and b a e not known. Instead they are 
approximated by 


K > 2 °b 2 


2 • /a ®\2 2 

°y = (Ax n } °b 


(2) Variances are carried instead of standard deviations. 

Thus, actual numerical values may be quite large. This 
means that the so-called round-off or over-flow problems 
arising from a limited word size computer need to be addressed. 

(3) Sometimes it is useful to limit the values of o x (o y ) between 

a minimum and a maximum. This will prevent the gains from 
becoming too high or too low. 
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Table 3 shows the Kalman gains and estimation error standard deviations 
depending on the input noise magnitude, o fl . The measurement error magnitude 
was 825 ft. When Table 3 is compared with Table 2, the general trend is 
clearly the same; i.e., gains become larger as the sampling time becomes 
longer. The proportional gain (o) seems to be very close (especially the 
o *32 fps case) ; however, the Kalman velocity gain is at least 2 to 5 
times smaller. This implies that using the Kalman filter algorithm to 
determine the filter gains is not unreasonable, if the acceleration noise 
magnitude was used to tune the performance. 

Figure 1 shows the "steady state" Kalman gains as functions of range. 
Along side the gains, the measurement and estimation error standard de- 
viations are plotted. Figure 1 (a and b) corresponds to the bearing error 
of 1 deg and 1 and A sec sampling periods. Figure 1 (c and d) corresponds 
to the error magnitude of 2 deg at the same sampling periods. As can be 
seen, the x error is proportional to range, whereas the corresponding esti- 
mation error magnitude becomes large at the far range, the filter gains are 
lowered proportionally to reduce the "total" error magnitude (balancing 
the acceleration uncertainty and the position uncertainty). 


Method (c) - Table-look-up Gains This method is a compromise between 
the last two methods. The first method is an optimization based on a 
numerical minimization procedure. The Kalman algorithm is more of an analy- 
tical method. The art resides in choosing the input (acceleration) uncer- 
tainty. Basically, the filter gain optimization depends on four parameters: 
(1) filter time constant, ; (2) sampling period, At ; (3) measurement 
error magnitude, o ; and (A) acceleration input uncertainty, a . That is, 
the estimation error is given by functions of the form 


x rms 

» 

x 

rms 


f x (T f , At, o x , ; x ) , 

f 2* T f» At ’ °x* a x^ 


(30) 


Thus, the objective is to minimize x ra8 (or * raB ) with respect to the filter 
time constant, Tf, from which gains may be computed. The difficulty is that 
f^( a ) and f 2 ( a ) ate not analytically tractable for a minimization procedure. 
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Table 3. Kalman Gains and Error Standard Deviations 


= 825 ft, » 8 fpss 


e/At 

-l 


°x (ft) 


0 ^ (fps) 


At (sec) 


6/At 

-1 


o (ft) 

X 


(fps) 


2 

3 

0.247 

0.272 

0.015 

0.014 

410 

430 

34.1 

33.5 


2 

3 

0.282 

0.362 

0.023 

0.027 

438 

496 

53.4 

55.4 


5 

6 

0.372 

0.414 

0.017 

0.018 

503 

531 

35.0 

35.6 



* 825 ft, a =16 fpss 

X 3 


5 

6 

0.482 

0.53 

0.031 

0.033 

573 

601 

57.6 

58.3 



= 825 ft, o =32 fpss 

X £L 


At (sec) 


6/At 

-1 


°x < ft > 


(fpt) 


2 

3 

0.374 

0.47 

0.043 

0.049 

505 

566 

88.3 

91.1 






































































































Figure 1. “ Continued. 



The gain optimization process can be facilitated further, if the gains 
a and 6 are expressed by a single parameter. For this reason, we confine 
this analysis to the critically damped gain configuration. With this 
assumption, the gains are given by 


a 

6 



(1 - Y) 2 • 


(31) 


Here, the parameter y is related to the filter bandwidth (■ I/t^) by 
the expression 

Y - exp ■ exp (-A T /T b ) . 


Furthermore, to account for the range effect, the position measurement 
error is assumed to be 


o 

X 


r °b • 


(32) 


If the expressions (31) and (32) are substituted into the rms error 
Eqs (25) , the following are obtained 


rms 


T A a „A 

Y At 
_(1-Y>' 


At 4 .2 . (1 -y) (1+Ay + 5y 2 ) _2_ 2 
, a + ^ r 


1/2 


(1+Y) 


7 

X 


rms 


[ll+Mi Ati ~2 + A (l-Y) 3 

L(l-3 ' 


-Y) 


(1+Y) 


1 2 2 1 
17"* J 


1/2 


(33) 


The above expressions can be optimized for Y(or t fe ) in terms of the sampling 
period and range if proper values of a and are given. The importance of 
At and y is that they are filter operating parameters, whereas a and are 
filter design parameters. 

Figures 2a through 2c show the plots of T fe * (quantized to a nearest 
second) in terms of and At which minimizes x^^. Figure 2a corresponds to the 
Non-aided configuration with a • 16 fpss; Fig 2b correspond? to Own- 
aided configuration with e ■ 8 fpss; and Fig 2c corresponds to Own-and 
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Target-aided configuration with a ■ 2 fpss. When the acceleration uncer- 
tainty is high, then short. When the measurement error is high, 

T b * is long. When the sampling time is longer, then so is t^. Depending 
on the filter configuration, can be selected and interpolated with a 

simple two dimensional table-look-up procedure. It may be desirable to 
"smooth" the time constant before computing the gain values. This will 
prevent any occurrence of abrupt change in the filter gains. 

It may turn out that this method may require as much computational 
load as the Kalman filter algorithm Eqs (24) through (27) when it is real- 
istically implemented. In the next section, preliminary simulation re- 
sults based on Monte Carlo passes are presented for fixed and Kalman gain 
filters. The results for Method (c) fall somewhere in the middle of 
these cases. 



-Aided Configuration 






Preliminary Monte Carlo Simulation Results 


A Monte Carlo simulation was set up to compare and analyze the per- 
formance of various estimation configurations. Because the impact of 
dynamics due to turn maneuvers is the significant key to filter performance, 
the various configurations are tested mostly during maneuver. Table 4 
lists pertinent kinematic parameters for Own and Target aircraft. Appendix 
B gives a detailed description of the aircraft dynamic model used through- 
out the simulation study. 

Table 4. Kinematic Parameters for Simulated Aircraft 


— — 

Own 

Target 

Ground Speed 
(kt) 


382 

Bank Angle 
(deg) 

15 

-15 

Turn Radius 
(nmi) 

2.2 

7.2 

Turn Rate 
(deg/sec) 

1.5 

-0.85 


Figure 3 shows the horizontal fzjectories in an earth fixed coordinate 
system. Figure 4 shows the relative position and velocity with respect to 
an Ownship north reference system, i.e., the coordinate system is not ro- 
tating as Own aircraft heading rotates. Figure 5 shows the time plot of 
Own and Target headings and ground speeds. Both aircraft turn simultaneously 
at t - 35 sec with bank commands of 15 deg (Own to the right. Target to the 
left). The turns last 120 sec. Because both aircraft are pulling 0.3 g, 
the relative acceleration could be as high as 0.6 g (20 fpss). 

Original Tracker Performance Figure 6 shows the statistical time 
plots of position and velocity errors of the TCAS 06 tracker at the samp- 
ling periods of 1, 4 and 8 seconds. The statistics are computed based on 
sixty (60) Monte Carlo passes and shown as the mean, mean plus one sigma and 
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(nmi) 


y (knt) 


Figure 4. Relative Position and Velocity of Target 
with respect to Own. 



Target 


200 knt 






mean minus one sigma, i.e., that approximately 60% of the errors would fall 
between the + sigma plots. 


Some observations are listed below. 


1) The range effect is clearly shown. As the target closes in range, 
the errors get smaller; as the target flies away, the errots 

get larger. 

2) The mean position errors remain relatively small (within 250 ft) . 

The dynamic delay effect of turn maneuvers does not seem to be 
significant. This implies that the position gains are on the 
high side. The maximum rms error is approximately 1000 ft at 

16 nmi range. Position errors tend to increase as the sampling 
time increases. (This was predicted by the theoretical results.) 

3) The velocity errors show the same range effect characteristics. 

The mean errors show the marked dynamic delay error especially 
for the y-axis. It is apparent from Fig. 6 that the y-axis 
contains most of the dynamics. The velocity gains are well tuned 
in the sense that the dynamic error and the high frequency error 
have similar magnitude. The maximum rms error reaches approximately 
100 kt. Except for the initial transients, the error curves are 
similar. It is notfed that the initial velocity estimate using 
the first two consecutive measurements becomes more accurate 
as the sampling period becomes longer because of the greater 
signal-to-noise ratio. 

Effect of Limltors on the Feedback Signal In many cases of navigation 
filter or feedback regulator control system design, a limitor is used on the 
feedback signal to prevent "unreasonable" error signals disturbing the 
system. In equation form, a limitor is given by 


Ax ^ " ■ 


Lx , if | A~x| £ L x 

sign (Ax) L x , if |Ax| > L x 


(34) 


That is, the feedback signal Ax may be replaced by the term L x in 

the filter Eq (8) depending on the magnitude. The y signal is generated 

In a similar manner. 
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Figure 7(a) shows the Monte Carlo simulation results with 
L * o and L * o 

xx y y 


(35) 


at the sampling period of 1 sec. While the x-axis behavior remains 
"reasonable", the y-axis clearly shows the unstable nature mostly due to 
not being able to track velocity. Figure 7(b) shows the case when two 
sigma values are used as the llmitor value at At of 5 sec. In this case, the 
unstable nature becomes apparent at At of 4 sec or longer (At of 5 sec 
is shown) . 


Other nonlinear devices of this nature were tested with similar 
results. For example, a soft llmitor defined by a quadratic function 
(rather than a straight line linear function) was also tested. This 
indicates that the performance of the current tracker can not be improved 
very much with these devices. 


Kalman Filter Algorithm Gains The Monte Carlo Simulation program 
was modified to include three filter configurations. These are Non-aided, 
Own-aided and Target-and Own-aided configurations. Instead of fixed gains, 
the filters utilize the Kalman filter algorithm for computing gains 
( Eqs (27) through (29)). Figures 8 and 9 show the statistical error time 
plots of three filters side by side for the sampling period of 1 and 4 
seconds. Figure 8 shows the case with the gains set at high values and 
Fig. 9 shows the case with the gains set at low values. 

The gain computations are done automatically. But, by adjusting the 
acceleration uncertainty value, a & , high or low gain can be selected. 

The input acceleration uncertainties are shown in the following Table 5. 


Table 5. Input Acceleration Uncertainties 

(o.. or o- in f/s^) 
x y 



Non-Aided 

Own-Aided 

Target- and 
Own-Aided 

High Gain 

32 

32 

32 

Low Gain 

20 

16 

8 
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Figure 7. Continued 
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Figure 8. High Kalman Filter Gains for Three Filter Configurations. 
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The rationale of using successively lower o^'s for the low filter gains is 
that the uncertainty should decrease as Ow.i and Target onboard signals 
are used to complement the relative position measurement. Note that for 
the high gain case, the acceleration uncertainty value is the same for all 
three configurations. Thus, the gains would be identical. This implies 
that the smoothing characteristics of alternating the measurement noise 
magnitude would be the same* 

Now referring to Fig. 8, the + one sigma envelopes are very similiar for 
all these configurations confirming the intuition expressed in the last 
paragraph. Even though the gains are computed to attenuate the range effect by 
lowering the gains at longer range. the effect is still very much prominent. 

The performance of the x-axis of the non-aided configuration looks similar to 
a fixed gain case. rhe y-axis shows similar traits; however, the mean 
position and velocity errors are larger for the Kalman filter case. This 
points out that x- and y-axis filter gains are computed separately. This 
means that the Kalman filter band-width for the y-axis is narrower than the 
constant gain filter. The mean y error for the 4 sec case shows smaller value. 


For this particular case of relative kinematics (a counter example 
will be shown later), the Own data and Target and Own data complementation 
successfully improve the tracking performance. This can be seen as smaller 
mean errors. For the case of Target and Own data complementation, the mean 
errors are remarkably close to zero even for the 4 second sampling time case 

Referring to Fig. 9, the following observations can be made : 

(i) Performance of the x-axis for both Non-aided and Own-aided 
configurations deteriorates. The major cause for this 

is the larger mean errors due to dynamic delay. Performance 
of the Own-aided configuration is larger than that for the Non- 
aided configuration. For the y-axis, the performance of both 
configurations are very similar. These comments are applicable 
for At ■ 4 sec case also. 

(ii) For the Target- and Own-aided configuration, errors are 
generally smaller. Also, the mean errors are small 
regardless of sampling time. The range effects are very 
slight in the velocity errors. This seems to indicate 
that the position measurements are used mainly to update 
the low frequency error in the average acceleration input 
and the high frequency input error is "integrated" out. 
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Non— Aided and Own-Aided Configurations it could be concluded that a 
partial indication of maneuvers is better than none. In particular, it seems 
true that complementing with Own data, even though lacking Target data, is 
better than no complementation at all. This would be true in general cases 
of the relative kinematics involving only Own aircraft maneuvering. In the 
case of both aircraft maneuvering, the filter performance depends on relative 
acceleration. For example, in the example used in the previous sections, the 
Own-aided filter performed better than the non-aided filter if the filter 
gains are the same. 

A counter example can be given wherein this is not the case in that 
instead of continuing to turn at the given bank angles, both aircraft go 
into bank reversal maneuvers simultaneously at the mid-point. Figures 10 
through 12 show the pertinent parameters. It is noted from the relative 
velocity plot of Fig. 11, that the relative x velocity is nearly constant. 
Therefore, in this case, the non-aided configuration should perform better. 
Figure 13 shows the error time plots of the Monte Carlo run. Both filters 
are set in the high gain Kalman filter mode. As can be seen, the position 
errors are comparable. The relative x— velocity estimates are markedly 
different : The non-aided filter shows a small mean error; whereas, the 
Own-aided filter shows a substantial non-zero mean error. The fact that 
the standard deviations are similar In size indicates that the differences 
in characteristics are due solely to the difference in (absence or presence 
of) acceleration input. 

Conclusions 

Important conclusions are summarized below: 

(1) Combination of Target-and Own-data complementation and the 
Kalman filter gain computation exhibited the best results. 

One advantage is that the errors due to the dynamic lag in- 
duced by maneuvers are non-existent; 

(2) Own-data complementation helps in cases where only Own is maneuver- 
ing. If both Own and target are maneuvering, then the estimation 
performance depends on the relative acceleration; 
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Figure 12. Time Plots 
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(3) Non-aided configuration with the constant gains has problems when 
either aircraft undergoes a maneuver. However, as will be shown 
later, its performance is credible when the relative kinematics 
is rectilinear. Because the Kalman filter gain computations do 
not require significant real time, this method rather than fixed 
gains is recommended; and 

(4) Adaptive feature of monitoring the innovations (or measurement 
residues) to detect dynamic delay error was not considered. This 
feature may be necessary to avoid using estimates which contain 
large and sustained dynamic lag errors. Without such a feature, 
the system does not know "when it does not know". 


The Impact of these errors is discussed in Chapter V with respect 
to CDTI and CAS applications. Also in the next chapter, their inpact is 
discussed with respect to the range and range rate estimate derivations. 
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III 

RANGE FILTERS 
Introduction 


By far the two most important kinematic variables for the TCAS system 
are relative range and altitude. These are monitored very carefully for 
targets in the vicinity of Own aircraft to determine separation status. The 
monitoring process involves essentially three steps: measurement , filtering 

and collision avoidance logic. The filtering algorithm provides the dynamic 
state estimates to the collision avoidance logic based on the measurements 
obtained by the surveillance subsystem. 


Note: The planned TCAS II software architecture incorporates similar 

filter algorithms in two different submodules, i.e., the filter 
algorithms are not separate from the other two functions. One set of 
filters is used within the surveillance subsystem. It is used to 
correlate and distinguish the internal track files and incoming 
measurements. The so-called "gating" techniques are used for this 
purpose. The gate's upper and lower thresholds are dynamically 
generated based on the predicted state variables. 

Another set of filtering algorithms are incorporated within 
the collision avoidance logic. Of course, these are used to provide 
the threat assessment. In our discussion, the filtering algorithms 
are treated as an independent functional block. 

For the range axis, the draft TCAS II Minimum Operational Performance 
Standards [9] proposes three filtering methods: 


(1) the horizontal x-y filter ; 

(2) the range oB tracker ; and 

(3) the range squared aBy tracker. 


These are discussed in the following sections. 


Range Filter Algorithas 



Range and Range Rate Estimates Derived from Horizontal Ax-Ay 
Estimate This method is recommended in the draft TCAS II MOPS when 
reasonably accurate relative bearing is obtainable. The algorithm is 
simple. The relative range Ar is given in terms of the Ax-Ay horizontal 
components by the formula 

Ar - I Ax 2 + Ay 2 ]*^ 2 . (3 

Taking time derivative 

Ar * lAx 2 + Ay 2 ] ^ 2 (Ax • Ax + Ay • Ay) 

* (Ax • Ax + Ay • Ay)/Ar . (3 

The range and range rate estimates can be obtained by substituting the 
horizontal estimates obtained in Chapter II. 


The range and range rate errors are related to the horizontal x-y 
errors (within the linear term) by 

m. • «. 

Ar » cosb Ax + sinb Ay , 


Ar 


X r • 

lsinb [sinb Ax - cosb 
+ cosb [cosb Ay - 


+ cosb Ax + sinb Ay 


• „ 

Ay] Ax 

sinb Ax]Ay) 


(3 


where 

cosb 


Ax 

Ar 


sinb 


A l 

Ar 


(3 


Equations (38) simply state that the range and range rate errors are 
the same order of magnitude as the respective horizontal components. 


The Range aB Tracker In this configuration, the measured range 
is fed into an aB tracker algorithm. The filtering algorithm is 
repeated below. 


range prediction: 


Ar + ■ Ar + At • Ar 
n n 

ID A + 

range error: Ar ■ At q+ ^ “ Ar 

(40) 

* + _ 

range update: ^ r n+l * + a 

range rate update: ^ r n+l * ^ r n + 

When the range measurement is missing, the predicted range is used for 
a short duration of time (6 sec) before the surveillance process is re- 
started again. The first two consecutive measurements are used to 
initialize the states. The recommended a and 8 gains are 0.67 and 
0.25 respectively for the nominal sampling interval of 1 sec. 


T he Range-Square aBy Tracker One of the major drawbacks of the range 
08 tracker is that the range is a nonlinear function of time even though the 
underlying relative kinematics is rectilinear. For example, if 

Ax * AXq + Ax • t ; Ay ■ Ay^ + Ay • t 


with constant Ax and Ay, then the range is given by 


Ar » [(Ax 0 2 + Ay Q 2 ) + 2(AXq Ax + Ay Q Ay)t + (Ax 2 + Ay 2 )t 2 ] 


1/2 


(41) 


The above expression can be approximated by the linear expression 
Ar * Ar^ + Ar t , 

which forms the basis for the aB tracker formulation, if Ar is reasonably 
constant. See Eq (40). 


Equation (41) forms the basis of the range square filter. 

both sides of Eq (41), one obtains a quadratic equation 

2 • 1-2 
s i Ar - s Q + s Q t + 2 Sgt » 


Dy squaring 


(42) 


51 


where 


8 o “ * x o 2 + * y o 2 » 

8 0 - 2 C&X Q Ax + Ay c Ay) , 

f A * 2 , A * 2 v 

s Q - (4x + Ay ) . 


• «> 

The quadratic coefficients s q and s^ will be constant as long as 
• • 

AXq, Ayg, Ax and Ay are constant. 


Because there are three unknowns, the underlying state equation must 
be three dimensional. The state equation is given by 


s 


1 At At 2 /2 


s 

• 

s 

- 

0 1 At 


• 

s 

s 

n+1 

1 

o 

0 

1 




(43) 


The observation equation can be obtained simply by squaring the range measure' 
ment, i.e.. 


s n+l 


< 4 '«/ 


(44) 


Equations (43) and (44) are in the form suitable for developing the so-called 
agy tracker algorithm. The algorithm is given below. 


prediction: 



1 At At 2 /2 


s 



* 

0 1 At 


S 

00 1 


A 
• • 

_S_ 


measurement error: s 


z . m .2 4- 


state update: 


A 


+* 1 



8 


s 


a 

9 

S 

- 

s+ 

+ 

B/At 

_*8^ 

n+l 

_‘s + _ 


_y/At 2 _ 


( 45 ) 


I 

i 

j The states are Initialized using the first three consecutive 

J measurements as follows 

m 

- s 2 . 

- 13 sj - 4 sj + sJ]/2At , <* 6 > 

* Is™ - 2 s™ + s™]/At^ 

8 are time scheduled according to 


* 

Sr 


The gains 


t 

? 


a 

n 


/ 3 (3n 2 + 3n + 2 ) 

) 

l * 


! 18(2n+l) 
D 

•I 



where 

D - (n+l)(n+2)(n+3) . 

n 


for n ■ 3 to 15, 
for n > 15. 


for n * 3 to 15, 
for n > 15. 


(47) 


for n ■ 3 to 15, 
for n > 15. 


Figure 14 shows the time schedule characteristics for the three gains. 

All three gains begin with large values and decrease to smaller values as 
more measurements are incorporated in the estimates. The gains are 
prevented from becoming too small so as not to lose the filter "adaptability . 
These are characteristics exhibited by a Kalman filter. 


As usual, when a measurement is missed because of surveillance failure, 
the predicted position, s + , is used in place of the measurement (without 
advancing the gain computation time frame). 
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where Ar is the true range and S r is presumably an independent random 
number. Then, the above becomes 

s® - Lx 2 + 2&r£ r + C r 2 “ s + C s . 

That is, the measured value is the sum of true value and 'measure- 
ment noise", which sort of satisfies the usual filtering assumption. 

However, the noise, £ g , is no longer independent. It depends on the 
state, s. Therefore, the proposed a$y tracker algorithm may not be 

optimal. 

Correlated Error aS Tracker Up to now, we have assumed that the 
range measurement error is white noise, i.e., independent from one 
sample to the next. However, there is strong evidence that the error is 
strongly correlated [10]. If the low frequency (or bias) error term is removed, 

the error is governed by the following linear equation 


5 r,n+l 


P 


>r »n 


+ 


n rH-l 


(49) 


Here, 

p - correlation parameter with the empirical value of 0.6801, and 
n * a zero mean stationary Gaussian noise with one sigma value 
of 69.5 ft. 


i 


The steady state variance of the range error is given by 
a 2 - (l-p 2 )” 1 o 2 - (94.5 ft ) 2 . 


(50) 


When the measurement error is white noise, it is true that the <*B 
tracker given by Eq (40) performs near optimum when the variation in the 
range rate is negligible. However, the above is no longer valid when the 
error is known to be correlated, as in our present problem. A method is 
available to modify the original algorithm in order to take advantage of 
the fact that the error is correlated. The method is based on the optimum 

filter derived by Tarn [11]. 
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The modification involves one additional term in the feedback signal 
Tba modified closed loop filter equation is given by 



where p is the correlation parameter of the noise process. 

It is noted that the modification requires one additional memory 
cell to store the previous range measurement. 

It was shown that the above formulation does indeed result in 
performance improvement 1 12]. In a simulation study conducted earlier, 
the range rate estimation error was reduced by 23%. Furthermore, the 
formulation is fairly robust with respect to a small variation in p. 



Monte Carlo Simulation Results 


Simulation Scenarios Four of the previously discussed range filter 
algorithms were implemented in the Monte Carlo Simulation Program to 
obtain statistical performance data for the range and range estimation 
errors. The four algorithms are: 

(1) Non-aided x-y fixed gain tracker; 

(2) Target- and Own-aided Kalman filter; 

(3) Range a 8 tracker; and 

(4) Range square oBy tracker. 

Configuration (1) uses the proposed fixed gains (a * 0.67 and 8 = 0.25 for 
At * 1 sec. Configuration (2) uses the Kalman gain update formula with the 
acceleration uncertainly of + 16 fpss (+1®)* Nominal measurement errors are 
assumed to be white noise processes with the standard deviation values of 75 
ft and 1 deg for range and bearing respectively. 

Four typical encounter scenarios were chosen motivated by the 
collision avoidance logic applications shown in Fig 15. These are (a) 
tail chase, (b) route crossing, (c) head-on, and (d) parallel turn-in en- 
counters. Figure 15 shows the horizontal projections of each geometry in a 
north-east coordinate system. Own flies due north with a 200 kt ground speed 
starting from the origin. The target kinematic parameters - ground speed 

(V_) , heading (40 and the miss distance (M.) - are listed in the following 
G a 

Table 6. 

Statistical data were obtained based on sixty passes for each of the 
above encounter cases with nominal measurement errors for the nominal 
sampling period of 1 sec. To test the sensitivity of the estimates to 
measurement errors, the route crossing and parallel turn-in encounter cases 
were repeated with twice the nominal measurement errors; the head-on 
encounter was repeated with twice and four times the nominal values. 

Except for the parallel turn-in case, these encounter scenarios are 
rather benign compared to the scenario used in Chapter II. However, these ar> 
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Figure 15. Relative Kinematics for Test Cases 
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Table 6. Target Kinematic Parameters 


Encounter 

V 

a 

(kt) 

(deg) 

M d 

(nmi) 

Turn Rate 
(deg/sec) 

Case No. 

Figure No. 

Tail Chase 

300 

30 

0 

0 

a 

15,16 

Route 

Crossing 

200 

270 

0 

0 

b.l 

b.2 

17,18 

19,20 

Head-on 

i 

160 

140 

0.3 

0 

c.l 
c.2 
c. 3 

21,22 

23,24 

25,26 

Parallel 

Tum-in 

200 

0 

0 

1.5 (target) 
-1. 5(0wn) 

d.l 

d.2 

27,28 

29,30 


thought to represent typical CAS encounters. The proposed parallel turn-in 
scenario in the TCAS MOPS seems to be very severe. The scenario requires 
both Own and Target to perform 65 deg bank angle maneuvering at 600 kt 
ground speed resulting in more than 2g lateral accelerations. Therefore, 
this scenario was scaled down to reflect more representative values of 
commercial operations. 

Simulation Results: Tail Chase Encounter Figure 16 shows sample 

time plots of range and range rate estimates together with the true 
values. Figure 17 shows the statistical time plots of error mean and 
mean +lo . Thus 68% of the error would fall between the dotted curves. 



* 




i 


Range and range rate estimates derived from the x and y components 
of both non-alded and aided filters show very similar characteristics. 

The range errors are 44 and 41 ft and the rate errors are 6.8 and 6.4 kt 
respectively. The results follow the fact that the filter gains are 
comparable for both filters. The range rate error pulses appearing at 
approximately 100 sec are caused by division by very small range estimates. 


i 

i 

> 
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The range ot& tracker shows somewhat larger errors -57 ft and 18.2 kt 
for the range and range rate. The rate error "pulse" at t=100 sec is not 
caused by division by sioall numbers. It represents dynamic delay error 
due to the sign change in the range rate at that particular time. 

The range square otBy tracker shows the steady state errors of 80 ft 
and 9.0 kt. When the range is zero, the range rate error shows a peculiar 
doublet behavior. The effect of the gain reduction schedule is very 
apparent in the initial transient of the rate estimate. The filter settles 
down to the steady state operation at t “ 20 sec as expected. 

Route Crossing Encounter Time plots of the simulation results, 
for the route crossing encounter, are shown in Figs 18 through 21. The pre- 
vious comments generally apply to this case also. The error statistics are 
very similar to the previous case. 

Figures 20 and 21 show the results with the measurement error magni- 
tudes twice those of the nominal case. Compared to the nominal case, the 
error standard deviations are twice as large. 

The non-aided and the aided configurations show similar mean error 
characteristics in range with a maximum of 60 ft. Since the other two 
configurations do not show a similar symptom, the mean error is caused by 
the y axis component. Because the target track is due west, whereas 
Own's is due north, the bearing error affects the y axis directly which 
induces the range effect. 

The singularity problem (of dividing by small range) becomes worse 
for the range square tracker as the noise level becomes larger. This 
implies that the range square filter is not reliable at extremely 
small range. The range aB tracker suffers the same reliability problem 
caused by an entirely different effect - the dynamic delay error build 
up. One other problem encountered for the range square otBy tracker 
is the .proolem of u e "range square" state becoming negative. This 
causes the computer implementation problem of taking the square-r ot of 
a negative number. 
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(c) Range a8 Tracker 


(d) Range Square oBy Tracker 


Figure 18. Sample Time Plots of Range and Range Rate Estimates for the 

Route Crossing Encounter. 




(c) Range aB Tracker (d) Range Square ciBy Tracker 

Figure 19. Statistical Time Plots of Range and Range Rate Errors for the 

Route Crossing Encounter. 
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Head-on Encounter Figures 22 through 27 show the results of the 
head-on encounter simulation results. The measurement errors were varied 
twice and four times the nominals. The largest level represents the lower 
performance enhanced TCAS II capability. As before, the previous comments 
apply here mostly. Because the encounter does not involve a collision (the miss 
distance is 0.3 nmi) , the singularity problem does not manifest in the 
estimates. The 06 tracker shows a considerable dynamic delay error when 
the range is closed to the mlninum. 

As can be expected, general performance of the estimation algorithms 
degrade as the measurement noise levels increase. In fact the standard 
deviations are remarkably proportional. The range errors varied from 
44.4 - 173.4 ft for the aided configuration to 54.7 - 228.1 ft for the 
aB tracker. The rate errors vary from 6.8 - 28.5 kt for the aided con- 
figuration to 18.5 - 73.2 kt for the aB tracker. The former provided 
the most accurate estimates; the latter provides the poorest estimates. 

The range effect is apparent in the x— y based estimates, but it is not 
in the other two configurations. This should be intuitively understood 
because the latter two do not contain the bearing measurement whereas the 
first two do. 

In these rectilinear trajectory cases, the estimates based on the 
non-aided and the Target- and Own-aided configurations as well as the range 
square aBy tracker do not show substantial qualitative difference (except 
the singularity problem of the last tracker). The aided estimates show 
the smallest errors. Thus, if the underlying kinematics are rectilinear, 
then all three configurations would be equally effective. This follows 
from the fact that the filter model equations are exact; hence, no in- 
duced dynamic error. 


67 


(c) Range aB Tracker 


(d) Range Square a0y Tracker 


Figure 22. Sample Time Plots of Range and Range Rate Estimates for the 

Head-on Encounter. 
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(c) Ranee aB Tracker 


(d) Range Square aBy Tracker 


Figure 25. Statistical Time Plots of Range and Range Rate Errors for the 

Head-on Encounter (twice the nominal errors). 
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(a) Non-Aided x-y aB Tracker 


(b) Target- and Own-Aided 
x-y Kalman Filter 


e 



(c) Range aB Tracker (d) Range Square aBy Tracker 

Figure 27. Statistical Time Plots of Range and Range Rate Errors for the 

Head-on Encounter (four times the nominal errors). 






Parallel Turn-In Encounter Figures 28 through 31 show the simulation 
results for the parallel turn-in encounter. The encounter scenario involves 
both aircraft turning with 15 deg bank angle maneuvers. The combined ac- 
celeration is 0.54 g. Figure 28 shows that the range rate is no longer con- 
stant , and the range is not rectilinear. 


Now, by referring to Figs 29 and 31, the following comments can be made 

(i) The non-aided x-y tracker shows substantial mean errors with 
the maximum range error of 91 ft and the range rate error of 35.1 kt. 

The standard deviations are 93 ft and 15.0 kt for the nominal measure- 
ment error level. Surprisingly, the error magnitudes do not increase 
proportionally. This may be due to the fact that the number of the 
Monte Carlo passes are too small to provide the statistical precision. 
Note that the error curves contain more "roughness" for the higher 
error noise. 

(ii) The aided configuration does not show outstanding mean errors. 

The standard deviations are comparable to other encounter cases. 

(iii) The mean errors for the tracker are substantially smaller. 

The standard deviations are comparable to other encounter cases. 

This implies that the gain values are selected more on the basis of 
"tracking" rather than on the basis of "smoothing". During the 
maneuver, this tracker performs better than the other two non-aided 
conf igurations . 

(iv) The mean errors of the range square tracker show interesting 
features. The maximum mean errors were 82 ft and 159 ft for the 
range and 29 kt for the range rate. The standard deviations do not 
show substantial increase. The first portion is caused by the dynamic 
delay not being able to catch up with the acceleration. The mean 
errors show oscillation. The second portion is caused by the effect 
of low gains, i.e., the tracker can not catch up because the low gains 
prevent rapid recovery. 
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(a) N on-Aided x-y aB Tracker (b) Target- and Own-Aided 

x-y Kalman Filter 



(c) Range otB Tracker (d) Range Square a By Tracker 


Figure 29. Statistical Time Plots of Range and Range Rate Errors for the 

Parallel Tum-in Encounter. 
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(c) Range aB Tracker 
Figure 30. Sampl 


(d) Range Square ciBy Tracker 

Time Plots of Range and Range Rate Estimates for the 
Parallel Tum-in Encounter (twice Nominal) . 
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Conclusions 


Table 7 shows the summary of the four range estimation algorithms 
for the various simulation scenarios. The first three cases involve 
recilinear encounters and the last one involves simultaneous parallel 
turn-in maneuvers by both aircraft. The following comments apply: 


(i) The aided configuration showed the best performance in 
all cases. The estimation errors are similar in magnitude 
independent of encounter geometry. The errors are proportional 
to the input error magnitudes; 

(ii) The range and range rate estimates derived from the non-aided 
x-y filters were better than the other two range axis filters, ex- 
cept the parallel turn-in encounter. The major problem is the 
large and sustained mean errors caused by the filter dynamic delay; 

(iii) The range square filter performed credibly. The low gain 
nature is apparent in the initial transient error behavior and the 
dynamic delay errors; and 

(iv) The range a& tracker suffers from the nonlinear range behavior at 
or near the minimum range, even for rectilinear encounter cases. How- 
ever, the transient periods are relatively short due to its high gain 
nature. On the other hand, this high gain nature passes a large 
portion of the high frequency noise achieving less smoothing than 
other filters. 

Effects of these errors will be analyzed with respect to CAS or 
CDTI applications in Chapter V. As mentioned earlier, the range and range 
rate estimates play a very important role in detecting threatening targets. 
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Table 7. Performance Summary Table of Various 
Range Estimation Algorithms 


Encounter 

Non-Aided 

Target- and 

Range 06 

Range Square 




Own-Aided 


o 6 y 




40.5 

56.6 

49.5 

Tail Chase 

6 . 8 1 (2) 





6.4 

18.2 

9.0 


lx 

42.3 

36.2 

56.0 

49.5 


6.7 

5.2 

18.4 

7.5 

Route 






Crossing 

2 x 

94.0 

84.0 

113.6 

94.2 


12.5 

9.6 

35.9 

14.1 


lx 

53.9 

44.4 

54.7 

48.4 


10.3 

6.8 

18.5 

7.6 

Head-on 

2 x 

l . ■ - - — 

107.7 

87.0 

113.2 

101.0 



20.3 

14.3 

36.2 

15.8 


4x 

194.6 

173.4 

228.1 

191.6 


37.0 

28.5 

73.2 

31.9 



62.9 (91.5) (3) 4 

39.7 

58.1 

57.8 (81.9) 

Parallel 

lx 

15.0 (35.1) (A) 

6.1 

21.1 (18.0) 

14.7 (29.2) 

Tum-in 

2 x 

99.8 (92.9) 

80.6 

117.2 

105.1 (109.2) 


19.8 (37.3) 

11.7 

38.1 (14.8) 

21.1 (29.2) 


( 1 ) steady state standard deviation of range error, in ft. 

( 2 ) steady state standard deviation of range rate error, in kt. 

( 3 ) maximum mean range error, in ft. 

(4) maximum mean range rate error, in kt. 




altitude filters 


Because the altitude is the primary axis of separation for the 
collision avoidance logic, monitoring of vertical state variables is an 
extremely important function of the l'CAS. Therefore, the estimation 
algorithm is the key element providing the altitude and altitude rate 
estimates. 

Target altitude is decoded from Mode C or Mode S transponder replies. 

The measured value is a binary integer with the least significant bit 
representing 100 ft. Figure 32 shows a schematic diagram of the altitude 
measurement process. According to the current TCAS specifications. Own 
altitude can be measured in two ways - before or after the Mode C encoding 
process. Thus, the TCAS processor can access either the transponder altitude 
data (with the 100 ft quantization) or the TCAS altitude data with much 
higher resolution. Depending on which of these is used, a proper estima- 
tion algorithm needs to be chosen. 

The proposed vertical tracker algorithm for TCAS i sage is based on 
the MIT developed Level Occupancy Time (LOT) tracker. It was developed 
for the active BCAS application in order to overcome the 100 ft quantiza- 
tion of the encoded altitude measurement [13, 14). 



Pilot Display 


Transponder 
Altitude Data 


TCAS Altitude 
Data 


Figure 32. Schematic Diagram of Altitude Measurement Process. 



The basic idea of the LOT tracker is to estimate the altitude rate 
indirectly by estimating the time duration (called the level occupancy 
time) in a particular quantization level. If the altitude rate is 
constant, then so is the time duration. Thus, the estimate of level 
occupancy time, T, is given by 

T ■ T + k (T - T) , 0 < k < 1, (52 

meas 

where 

T ■ t - t • (53 

meas j ump last jump 

Then, the altitude rate estimate is given by 

* - (54 

z - 100/T 

Two of the ramifications of the LOT tracker algorithm are: 

(1) It requires at least two level changes to obtain rate infor- 
mat ion ; and 

(2) It requires at least three level changes to ascertain a 
change in rate (acceleration) . 

The vertical tracker software specification contained in the draft 
TCAS II MOPS is very complex owing to many heuristic logic elements. 

In the following sections two algorithms are discussed. These are 
( 1 ) two-state 06 tracker, and ( 2 ) level switching time algoritlm. 


Alpha Beta Tracker Algorithm 

Algorithm The aS tracker is a linear two-state recursive filter 
that estimates the aircraft's altitude and its rate of change based upon 
noise contaminated measurements of altitude. The equations for this 
algorithm are 





+ 

prediction: z * 

z + At • z 
n n 

altitude error: z m 

m + 

Z n+1 

altitude update: 

« z + + az 

altitude rate • 

update: n+1 

* z + (B/At)z 
n 

m 

where the measurement, z , 

is given by 

z m m z + ^ . 



( 55 ) 


(56) 


As usual, when the measurement is missing or invalid because of the 
surveillance failure, the predicted value is used in place of the 
measurement. The first two consecutive valid measurements can be used 
to initialize the estimate as follows: 


m 


1 ’ 


and 


<•; - z o )/At 


(57) 


Gain Selection There are many methods to select the feedback 
gains. A few methods were discussed in previous chapters. The basic idea 
is to tune the performance by compromising between the conflicting 
requirements of a fast response filter and of a good noise smoothing 
filter. A fast response filter would have a short time constant or 
vide bandwidth, whereas a good noise smoothing filter would be a sluggish 
system with long time constant or narrow bandwidth. Thus, it would be 
natural to end up with a different set of gains depending on different 
performance measures. For example, by optimizing the error due to step 
drift" in velocity, Benedict and Bordner 0.5] obtained a relationship 

between a and B , 


Now, a can be varied to satisfy other filter performance specifications. 
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The other method is to use the so-called exponential weighted least 
square error [16 ,17 ] . With this criteria, the gain values can be expressed 

. (59) 

or 

8 - 2 -a -2 /l-a . (60) 

Parameter y is related to the filter bandwidth by 

y * exp(-o>j j At) . (6l) 

The gain values recommended by the TCAS MOPS are 0.5 and 0.15 for 
a and 8 respectively, at the nominal sampling rate of 1 sec. 

These gain values are applicable when the altitude measurement error 
can be modeled by an independent random sequence. When it is no longer 
independent, as pointed out by the Billmann's study [10], The algorithm 
itself needs to be modified along the line suggested by the range filter 
with correlated measurement error. (However, the gain selection problem 
is a minor one compared to the problem caused by the 100 ft Mode C altitude 
quantization.) 

Figure 33 shows sample time plots of altitude rate istimates using two 
different sets of gains. The results were obtained by the a8 tracker with 
altitude measurement containing correlated additive error only. The 
measurement sequence was obtained by a level-climb (10 fps) -level flight vertical 
profile. (The experimental set-up will be explained in more detail later.) 


parametrically as 

. 2 
a • 1 - Y 

6 - (1-Y) 2 


84 




ifc 


5 ) 40 60 80 100 120 140 160 180 200 220 240 260 280 300 

Time (sec) 


i) ex = 0.5 and g = 0.15 Case 


o o o o Estimated 



Time (sec) 


(b) a " 0.486 and g = 0.0803 Case 


Figure 33. Altitude Rate Estimates for Different Set of Gains 
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These examples are used only to illustrate the effect of filter gains. The 
top plot shows the case with the recommended gains of a * 0.5 and R = 0.15. 
The bottom plot shows the rate estimate using Eq (59) with y * 0.717. For 
the sample case, the latter is clearly better. 

Effect of Quantized Measurements The Traffic Alert and Collision 
Avoidance System requires estimates of altitude and its rate of change 
for an intruder aircraft, based upon the Mode C altitude reports. These 
reports provide the altitude measurements to the nearest 100 ft. If a 
classical aB tracker is used, then a 100 ft "stair case" input induces an 
undesirable transient response each time the quantization level is crossed. 


This can be explained by noting that within an altitude level, the 
alLitude measurement, remains unchanged. This causes the estimate of 
vertical speed to decrease as estimated and measured values of altitude 
become the same. Figure 34 shows the estimates obtained by the aB tracker 
with Mode C reports as input. It clearly shows the undesirable train 
of transients. It is noted that the simulation result was obtained using 
a = 0.28 and 8 = 0.06. These values are used in the altitude tracker 
within the surveillance module rather than the CAS logic. 

The basic problem with the quantized measurements is the nonlinearity. 
The apparent error magnitude depends on the actual altitude as well as 
other additive error sources. The following two examples identify problems 
caused by quantization. 


Example 1: A 100 ft altitude jump occurs when the actual altitude 

goes from 10.049.9 ft to 10,050.1 ft. This represents an actual 
rate of + 0.2 ft/sec at the sampling rate of 1 sec. The- tracker 
will process this jump by increasing the rate estimate by 8 x 100 
ft/sec (6 is the rate gain). With a typical value of 0.15 for B, 
this implies that the rate estimate jumps by 15 ft/sec. 

Example 2: When an aircraft is flying at 10,050 ft altitude, then 

a small magnitude (say 0.2 ft rms) high frequency (random) noise 
will result in a situation where consecutive measurements are 
100 ft apart. 
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The above mentioned problems indicate that the capability of the aB 
tracker in response to quantized measurements is very limited. Alterna- 
tive algorithms must be considered, which is the subject of the next 
section. 

When the altitude measurement is obtainable with a much finer resolu- 
tion than 100 ft, then the aB tracker algorithm would be near optimum*. 

This is the case with the Own altitude input. Thus, the aB tracker (with 
a fine resolution altitude) can be used to judge other filtering schemes in 
terms of performance. 


* In many vertical filter implementations for navigation, guidance and 

flight control applications, other signals are available. These include 
the attitude angles, body rates and accelerations as well as the barometric 
altitude. Furthermore, these signals are available for processing at a 
much higher frequencey (50 msec compared to 1 sec for TCAS) . Therefore, 
the TCAS vertical aB tracker would not be able to perform as well as other 
implementat ions . 


True Velocity 




Level Switching Time Filter 


The unsatisfactory response of the aft tracker to the quantized 
measurements motivates the search for an alternative algorithm. As 
mentioned previously, the Level Occupancy Time (LOT) tracker algorithm 
was designed to fill this void. However, it has very complex software 
based on heuristic logic. The aim here is to reformulate the problem 
in a different light. 

The basic idea of the level switching time (LST) filter can be summarized 
in the following steps: 

1. Detection of level switching. 

2. Estimation of altitude and level switching time. 

3. Estimation of altitude and rate based on the last four 
estimates of altitude and level switching time. 

4. Modification of the estimated altitude and rate by error 
feedback. 

3. Validity verification of the modified altitude and rate values. 

The following sections explain each of the above steps. 

Leve l Switching Detection Logic The combination of measurement, 
noise and the quantization process can result in an erroneous 
indication of level change. To prevent such falsely reported level change, 
the decision on the level switching is made if three out of the past five 
measurements indicate such a change. This implies that with a sampling 
time of A, it takes at least 3A to detect a level change. This delay 
is necessary to compensate for erroneous level changes caused by noise 
and the quantization process, and any reduction in this time results in a 
reliability reduction of the level switching decision. 

Estimation of Altitude and Level Switching Time It is clear that 
more altitude information can be extracted from a segment of reports 
containing different levels rather than a single level. Consider, for 
example, a history of Mode C reports 
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Time 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

(sec) 

Report 

100 

100 

100 

100 

100 

100 

110 

110 

110 

110 

(100 ft). 


The segment corresponding to t»56 to 60 sec does not contain any informa- 
tion by itself other than altitude is between 9,050 to 10,050 ft. On 
the other hand the segment corresponding to t-59 to 63 sec contains a lot 
more information, viz., the true altitude was 10,050 ft some time during 
that interval and if the additive noise is small, then that time point would 
be between t of 60 and 61 sec. Thus, the latter segment is infinitely more 
useful in pinpointing the true altitude and the point of crossing that altitude. 
That is, the combination of true altitude and the corresponding time 
is the important factor. 


Motivated by the above argument, the following algorithm is used to 
estimate altitude and time: 




( 62 ) 


where z^ is the Mode C altitude report times 100 ft at t.. The four (4) 

1 l L 1 

most recent pairs of z and t are stored in shift registers. These are used 

to derive a rough estimate of altitude rate discussed next. See Fig. 35 

for illustration. 


Also stored is the estimated level, L, which is the quantized value 

, L . 
of z , i.e., 

L - Int[(z L + 50)/100] . (63) 

If the levels show a cyclic behavior (i.e., the latest level and the one 
before the previous level are identical), then level flight is declared. 

Altitude and Altitude Rate Estimation The stored level switching 
variables {z£, t£) : k«l,..,4) (with k*l being the latest) are used to 
generate rough estimates of altitude and altitude rate. The estimates 
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Figure 35. Illustrat 


are obtained using the principal of least squares, assuming that the altitude 
rate is constant during the time interval of interest. There are four 
cases to be considered depending on the number of stored measurements. 


(i) Altitude Initialization: on i y on e level switching is 

observed, then the estimates are initialized as 

z » z^ and z » 0 . (64) 

(ii) Altitude Rate Initialization: When two level switchings 

have been observed, then the altitude and rate can be initialized as 

z * z^ and z » (z^ - zlp/Ctj' - t^) . (65) 


(iii) Three Point Least Squares Fit: When three level switchings 

have been observed, then the altitude and rate can be estimated as 
follows. Because of the constant speed assumption, three linear 
equations can be written at t*t^ . 


z^ ■ z + (t^ - t^)z + ^2 * 
z^ - z + (tj - t^)z + z 3 , 


( 66 ) 


uh »m £ are errors to be minimized. Using the standard least 
squares method, the estimates at t ■ t^ are given by 



1 

A 2 + A 3 " (A 2 + A 3 ) 



L L" 

z 2 + z 3 

1 

N •> 

1 

m — 

D 

-(A 2 + A 3 ) 3 


v" 

+ S 3 Z 3 


where 



and 

D ■ 2(A^ - + Aj) 


(67) 
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(iv) Four Point Least Squares Fit: When four level switchings have 
been observed, the previous method can >-e extended easily. The 
estimates at t«t^ are given by 


z 

1 

A 2 + A 2 + A 2 -(A 2 + A 3 + A 4 ) 


*i + 

L L , L 

*2 + *3 + *4 

3 

Z 

D 

-(a 2 + a 3 + a a ) a 


V2 

+ A 3 *^ + A 4 z \ 


( 68 ) 

where 

, i - 2, 3, A, 

and 

D - 4(A 2 + A 2 + A*) - (A 2 + A 3 + A 4 ) 2 . 

It is noted that A^'s are "measured" level occupancy times. 


a 

The estimates z and z are computed on the basis of the most recent 
level switching time, tj. The time epoch t^ is inevitably at least 2-3 sec 
behind the current time, t. Thus, the current time estimates need to be 
extrapolated according to 

z(t) - z + (t-t^)z (6-) 

The altitude rate estimate remains unchanged. 

The least squares solution given by Eq (67) or (68) could be recast 
in a recursive algorithm. However, with the current computer technology, 

Eq (67) or (68) are straight forward. The recursive form may be advantageous 
for the loust squares fit using six or more data points. 

Fine Tuning by Error Feedback The previous methods apply only when 
a level switching is detected by the "3 out of 5" rule. At the time of 
detection, a new level switching time, t^, and the corresponding altitude, 
z^, are available. This time t^ could be very much behind the current time. 
Thus, there could be a big time gap between the estimate updates. 
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The time interval during which the Mode C altitude reports Indicate the 
same flight level provide "passive" information. The information will not 
allow us to compute the estimates, but it provides a means to check if 
the available estimates are still reasonable. This can be achieved by 
comparing the "predicted" altitude at the current time (given by Eq (69)) 
with the current Mode C altitude. If the difference between the two is 
less than 50 ft (a half of Mode C quantization) then the last estimates 
seem to be still accurate. On the other hand, a difference of more than 
50 ft Indicates (a) the last estimates are in error; (b) the target 
dynamics has changed; or (c) the Mode C report is in error due to noise. 

The estimates need to be updated in cases (a) and (b) but not in case (c) . 
This suggests that the average difference over some time period should be 
used rather than instantaneous difference. 


To accomplish this, the average value of error for the past five 
sampling times is calculated. Any amount in excess of + 50 ft is fedback 
to modify z and z. Therefore, z and z are updated according to the 
following equations: 


where 


and 


^ Z ^New 

(•'old + V 1 

/v 

• 

(z) New 

(z) 0ld + k 2 Z 

1 * 

- 50 , 

z “ ) _ 


t - 

+ 50 , 

- 1 

V' / c ~ \ 

Z ’ 5 

S ( z ‘ ■ J 


if z > 50 

if z < -50 


(70) 


(71) 


(72) 


The z^'s are the Mode C reported altitudes (times 100 ft) at the five most 
recent valid reporting times, t^. (The reports and times are stored in 
five tier shift registers and used in the level switching detection logic.) 

The z. are predicted altitudes extrapolated from the last level switching 
l 

time, t^ . Thus, 


‘Old + (t i " t i ) 


Old 


( 73 ) 
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a nd kj are the feed beck gain* for the altitude and rate update 
The valuea of and k 2 mu at be selected by compromising Detween the re- 
quirements for rapid modification and possible instability (divergence of 
i^’#) caused by large values. Values for k 2 of 0.632 and k 2 of 0.155 are 
reasonable values as vill be shown in a later section. 

Validity of the Modified Altitude and Rate Values To verify vali- 
dity, the modified values of z and * are used to recompute the estimates of 
altitude at the last five sampling times, and the average error z is re- 
computed. If the absolute value of the resulting average error is less than 
50 ft, the validity of the modified values is established. To avoid ex- 
cessive modification of altitude, this modification is limited to + 35 ft, i.e., 
the modified estimate z Ngw in Eq (70) would not be changed by more than 35 ft 
from the originally computed £ Qld by the least squares algorithm. 

Figure 38 shows an over-all computational flow of the proposed Level 
Switching Time filter. 


Remark A basic ingredient of the proposed tracker algorithm is 
the idea of obtaining the rough altitude and altitude rate estimates 

based on the estimated altitude and level switching time (z^, t ) pairs. 
The three-out-of-five rule was specifically designed for the sampling 
period of 1 sec. The same rule may not apply or be desirable for 
other sampling periods. Therefore, other methods of determining the 

(z , t^) pair must be devised. However, the rest of the algorithm 
should be applicable without modifications. The modular construc- 
tion of the logic would ease the task of replacing the level switching 
time detection logic. 
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Figure 36. Level Switching Time Filter Macro Flow Chart. 
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Monte Carlo Simulation Results 


This section describes the simulation set up, and it analyzes the proposed 
Level switching Time (LST) filter algorithm as well as the ag tracker. 

The ag tracker algorithm with non-quantized altitude input will be the basis 
of comparison. 


Simulation Set Up The aircraft vertical dynamics are simulated by 
a simple second order altitude select/hold logic. In equation form 



(74) 


where the acceleration command, z £ , is given by 


U — ( i c - z + tj)!) 

v ref 


Also, the velocity command, z £ , is given by 


z “ [[ — (z f - z) ]] • 

c t ref z , 

z ref 


(75) 


(76) 


In the above expressions the term £• is the altitude rate noise injected to 

z 

simulate such phenomena as wind gusts and pilot activities. It is designed to 
change magnitude on a 20 sec average basis, with a standard deviation of 
0.75 fps. The notation [[ x]] means x is authority-limited to y. z f , 
z re £ and z Te f are the altitude select references, i.e., they are the 
desired altitude, maximum rate magnitude and maximum allowable acceleration. 

a 

Once z re f is selected, the vertical velocity command, z c , is generated by 
Eq (76). If the current altitude is not close to the reference altitude, 

a 

the co nmand takes on the maximum value, z re f If Che altitude is close 
to the reference, then the command takes a smaller value. The acceleration 
command, z c , is generated in a similar manner depending on the (previously 

a a 

computed) velocity command, z c , and the current velocity state value, z. 
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1*9 + 


For example, if Che current velocity is close to the velocity command 
within the pilot noise input, then the acceleration command takes on a 
small value. Therefore, by selecting values for z rgf and * re{ (z ref is 
fixed to 2.5 fpss), reasonable altitude and altitude rate time histories 
can be generated. In the simulation program, Eqs (74)-(76) are computed 
at an integration cycle of 0.2 sec using a trapezoidal integration routine. 

The altitude generated by the above method is sampled at a one second 
interval. High frequency additive noise is added to the true altitude 
which in turn is quantized to the nearest 100 ft. Figure 37 shows the 
overall block diagram of this process. 

The additive error model is generated by a second order correlation 
process* as suggested by Billmann [ 10] . The generating equation is given 
by 

z * 1.066 i 1 - 0.191 2 _ + f . <77) 

n n-1 n-2 z,n 1 


Here, £ z is a random noise with a normal density function, mean value of 
zero, and standard deviation of 10.5 ft. By solving the associated 
Yule-Walker equation [ 18], it can be shown that the steady state standard 
deviation of the additive error is 23.8 ft. The measured altitude is given 
by 


m 

7 . 

n 


z 

n 


+ z 

n 


(78) 


represents the baro-al time ter output, and it is input to the ag tracker 
for generating Own altitude estimates. It is also input to generate the 
Node C reported altitude. 

The Mode C quantized altitude, z^ , is given by the equation 


* At a preliminary simulation stage, it was found that the level 
switching detection logic was sensitive to the nature of additive 
noise. For example, a two-out-of -three rule worked reliably with white 
noise error * it did not work for correlated noise. The three- 
out-of-five t works for both types of error. 
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Figure 37. Aircraft Vertical Dynamic and Measurement 
Process Model 
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[ z + 2 + 0.5 q 1 

JL—5 | (79) 

where q is 100 ft and Intlx] means the integer value of x. 

Figure 38 shows time plots of true altitude rate, altitude and the 
Mode C altitude reports (plotted every 3 sec) for the case of a pull-up 
maneuver of 2000 to 3700 ft at a rate of 10 fps. It is noted that the 
rate shows a small magnitude execursion about the command value of 10 fps. 

The aB tracker and LST tracker algorithms, together with the measure- 
ment generation module, were assembled together in a sixty pass Monte 
Carlo simulation program to generate statistical performance data. The 
true altitude and altitude rate time histories are computed once in the 
first pass. Afterwards, a new set of altitude errors are added to the 
true altitude which is fed into the tracker algorithms. 

Simulation Scenarios Several different simulation scenarios were 
examined during the course of this study. In all thirteen cases were 
run with the altitude profiles generated for vertical rates of + 5, 10, 20, 
and + 60 fps. Most of these are run to obtain raw statistical performance 
data. Some of these are run in conjunction with others to obtain "sensitivity 
data. Therefore, the total number of computer runs made during the study 
was considerable including the algorithm development, design and tuning stages. 


The following cases are discussed: 

1. Effect of quantization on the aB tracker; 

2. Effect of a and 6 gain variations; 

3. Performance change due to difference in additive errors; and 

4. Selected individual cases. 


100 



TIME (SEC) 


Figure 38. Time Plots of Vertical Rate, Altitude and 
Mode C Reported Altitude. 
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Effect of Altitude Quantization In the aB Tracker Algorithm Besides 
its use on the high resolution Own altitude input , the aB tracker may be 
used as a preprocessor in the surveillance function of correlating targets 
with the gating technique discussed previously. The algorithm i.e used 
with the Mode C reported altitude. Therefore, it is important to obtain a 
basic idea of how it performs under such circumstances. 

Quantization values of 25 ft and the regular Mode C 100 ft were 
tested. The quantization of 25 ft represents two more bits of altitude 
information in the Mode C transponder reply. Figure 39 shows sample 
altitude and altitude rate time plots for a 5 fps climb profile. The 
estimates are plotted every 3 seconds. Clearly, the quantization size 
has a large effect. For the 100 ft quantization case, the altitude 
estimate shows the stair-step characteristics. The rate estimate shows 
a sequence of transients, the est i mate jumps to large values and then decays 
to zero and so forth. For the 25 ft quantization case, the altitude esti- 
mate does not show the stair-step characteristics. It follows the true 
altitude fairly smoothly. The rate estimate tends to follow the true 
value more faithfully. The transient effect is still very much in evidence. 
This can be seen when the aircraft levels-off. In this phase, the rate 
estimate is jumping around. This is caused by the tendency for the measure- 
ments to change quantization levels more frequently. This indicates 
that even though the estimates are improved when the quantization is 25 ft, 
the classical aB tracker algorithm may not provide sufficiently accurate 
estimates. The following Table 8 summarizes the error statistics in terms 
of average rms. (Mean and standard deviation values do not mean very 
much when the estimates contain so much transient behavior.) 


Table 8. Quantization Level Effects on Error Magnitudes 


r:— - — — :■ — 1 ~ r- : - r, — i 

Quantization (ft) 

100 

25 

Altitude error rms (ft) 

35 

25 

Altitude rate error rms (fps) 

5.3 

4.3 

l -—Tl 
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Effect of the q0 Filter Gains The effect of filter gains is 
examined by comparing high and low gain configurations. The previously 
described correlated error model and a 10 fps climb altitude profile were 
used for simulation comparison. The high gain configuration used the 
o and 0 values of 0.5 and 0.15 respectively. This set is the recommended 
values for obtaining Own altitude and altitude rate estimates within the 
CAS logic. The low gain configuration uses the a and B values of 0.486 and 
0.0803, respectively. These gains are derived by the exponential weighting 
method. Sec Eq (59). The altitude measurements are not quantized. 

Figure 40 shows sample altitude and altitude rate estimate time 
plots. The altitude estimates are very similar, ref lecting the fact that 
the proportional (a) gains are practically equal. The altitude rate 
estimates are substantially different. The high gain configuration shows 
much larger errors compared to the low gain configuration. 

The reason for using a high rate gain (8) is to tr°ck the rate 
during the acceleration and deceleration phases, i.e., it provides a wider band- 
width. From the simulation results this point is not apparent. The 
low gain configuration does not show a particularly sluggish response 
when the altitude levels off. The mean peak error during this period is 
a little larger for the latter configuration, however. 

Table 9 summarizes the steady state error statistics in terms of 
error standard deviations. 

Table 9. Gain Effects on Error Magnitudes 










Tracker (b) Low Gain aB Tracker 








Different Error Characteristics The purpose of this section is 
two-fold. One is to exanine the effect of additive noise which is injected 
to measurement prior to the quantization process. The other is to com- 
pare the performance of the 06 tracker (without quantization) and the 
LOS tracker (with the standard 100 ft quantization). 

The error variations are (1) zero mean white noise with the standard 
deviation of 23.8 ft, (2) the correlated noise identified by Billmann 
(Eq (77)), and (3) the same noise but with 95Z report validity rate. The 
standard deviation of 23.8 ft for the white noise was chosen to match the 
steady state standard deviation of the correlated noise. The report 
validity rate is defined to be the probability of receiving a valid reply 
to an interrogation transmission. Thus, it means that 5% of the time, the 
Mode C reported altitude is missing. The 95X report validity rate is 
achievable with the current TCAS II design [19 ] • 

The flight scenario is a 20 fps climb profile from 2000 (FL20) to 4700 ft 
(FL47). Figures 41 and 42 show sample and error statistics time plots 
respectively for the white noise case. The aB tracker results are on 
the left and the LOS tracker results are on the right. Figures 43 and 44 
show the similar quantities for the correlated noise case and Figs. 45 
and 46 show results by using correlated noise and 95X report reliability. 

Table 10 shows the summary of the average statistics in comparing these cases. 

The following comments and remarks are derived from the simulation 
results. 

(i) As previously stated, the high frequency estimation errors 

(as shown by the standard deviations) are smallest for the white 
noise additive error case. For the aB tracker, the explana- 
tion is that this filter is based on the white noise assumption 
but not on the time correlation model. A most likely explanation 
for the LST tracker is that the level switching detection logic 
utilizes the sur of five measurements. In this case, indepen- 
dent errors tend to cancel each other. 
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Figure 41. White Noise Additive Error - Sample Time Plots 








Figure 42. White Noise Additive Error - Statistical Error Plots 








Figure 44. Correlated Noise Additive Error - Statistical Error Plots 




Figure 46. Correlated Noise Additive Error with 95% Reliability - Statistical Error Plots 



Table 10. Statistical Summary of Additive Noise 

Error Effect 



aB Tracker 

LST Tracker 


15.5 

(13.8)* 

16.7 

(85.7) 

Wh. e Noise 

2.5 

(12.9) 

1.6 

(20.8) 


28.9 

(23.5) 

31.2 

(83.9) 

Correlated Noise 

3.0 

(11.9) 

2.4 

(20.8) 

Correlated Noise 

28.7 

(18.5) 

34.9 

(85.4) 

with 




(20.5) 

55% Reliability 

4.6 

(12.3) 

3.3 


* Altitude standard deviation (ft) Peak mean error (ft) 

Altitude rate standard deviation (FPS) Peak mean error (fps) 


(ii) The mean errors for the LST tracker are very similar for the 
three cases. This indicates that the noise input has a 
secondary effect compared to the primary effect of the 100 ft 
resolution. The mean errors for the ag tracker are similar 
except for the altitude error for the white noise case which 
is smaller. 


(iii) The steady state altitude rate errors for the LST tracker 

are better than those of the aB tracker. The average improve- 
ment is 28%. During level light, the LST tracker is 
locked on to the nominal 0 fps. All of the LST tracker errors 
in this region reflect the rate noise injected to the altitude 
profile generation model. The altitude rate errors for the aB 
tracker are affected by the high frequency measurement noise. 

(iv) The aB tracker attenuates the white noise magnitude to 66% 
in the position estimate (15.5 vs. 23.8 ft). This is not 
the case with the correlated noise. In fact, the noise 
magnitude is amplified by 21% (28.8 vs. 23.8 ft). This 
implies that the filter is not tuned for the colored noise. 


(v) The high frequency altitude errors for the LST tracker are 

larger by only 12% compared to the aB tracker. This is despite 
the 100 ft quantization used for the LST algorithm. 
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(vi) The average peak errors for the aB tracker are 18.6 ft and 
12.4 fps for the altitude and rate. The coaparable numbers 
for the LST tracker are 85 ft and 20.7 fps. Essentially 
the LST tracker errors are. Induced by pure tine delay cor- 
responding roughly to 100/s. 


(vii)The dynamic delay times were 12.9 sec for the a6 tracker 
and 15.7 sec for the LST tracker. Here the dynastic delay 
time is measured by the average time interval between the 4 fps 
error point in the mean altitude rate error history. The delay 
is thought to be equivalent to a rise time of between 12.5 and 
87. 5X of the step response. 


(viii)The LST tracker altitude rate transient errors at the 
time of rate change are considered to be due solely to 
the quantization effect. For practical purposes, the same 
behavior occurs every Monte Carlo pass. This is evidenced 
by the lack of standard deviation spread during that period. 


(ix) When the altitude report reliability drops to 95% (i.e., 

one in twenty reports is missing, on the average) the estimation 
performance tends to degrade slightly. However, for 
practical intent and purpose, the degradation is not per- 
ceptable. 


Selected Individual Cases Three selected individual cases are 
presented and discussed in this section. The cases presented represent 
performances corresponding to different altitude rate profiles of -5, 

10 and 60 fps. (The 20 fps altitude rate profile was discussed in the 
previous section in detail. ) 

The extreme cases present particular challenges for the estimation 
task. The challenge for the 5 fps case is that the dynamic effect is 
lost or masked by the 100 ft quantization. New and useful information 
is available every 20 sec on the average. That is, the filter algorithm 
is "dead-reckoning" for 20 sec before the estimates can be updated. For 
the other extreme case of 60 fps, the challenge is that there is paradoxi- 
cally too much information. The level change occurs every 1.67 sec. Ther 
fore, small errors in the level switching time computation represents a 
large percentage of 1.67 sec; this induces errors in the altitude rate. 


Figure 47 and 48 show the sample and statistical error time plots 
for the -5 fps descent profile. Initial altitude is 2,970 ft, and final 
altitude is 2,030 ft. These intermediate altitudes were chosen so that the 
Mode C reports will be oscillatory in nature between, for example, FL2G and 21 
during the final level-off. The following comments and remarks are derived from 
the simulation results. 

(i) As seen in Fig. 47, the altitude rate estimate of the a$ 
tracker is very noisy, whereas that of the LST tracker is 
much more solid. The altitude estimates show a very similar 
characteristic. 

(ii) The initial time delay effect (-30 sec) of the LST tracker 

is clearly shown in the altitude and altitude rate estimates. 

This is caused by the basic 20 sec time delay until new and 
useful information becomes available. The LST tracker 
transient behavior at the final level— off is compounded by 
the intermediate altitude of 2030 ft. The Mode C reported 
altitudes are oscilatory, as can be seen in the altitude 
estimate. 

(iii) The standard deviations of rate errors are 1.8 and 1.4 fps 
for the aB and LST trackers, respectively. The peak mean 
errors are 3.9 and 4.9 fps. 

(iv) Time delay for the LST is approximately 30 sec. The transient 
error for the aB tracker is masked by the high frequency error. 


Figures 49 and 50 show the sample and statistical error time plots for 
the 10 fps climb altitude profile. Initial altitude is 2030 ft, and final 
altitude is 3730 ft. The following comments and remarks apply to this case. 


(v) Similar to the previous case, the LST tracker's estimates 
are more "solid" than those of the aB tracker. 

(vi) The time delay effect (-15 sec) of the LST tracker is clearly 
shown in the sample time plot of the altitude rate estimate. 
This causes an altitude overshoot (-100 ft) at the level-off 
period. 

(vii) The standard deviations of the rate errors were 2.0 and 1.3 
fps for the aB and LST trackers, respectively. The peak rate 
errors were 4.4 and 6.3 fps. The comparable altitude errors 
were 27.0 and 9.8 ft for the oB tracker and 25.8 and 52.? ft 
for the LST tracker. 
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Figure 47. Sample Time Plo. s for -5 fps Profile 
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Figure 50. Statistical Error Time Plots for 10 fps Profile 





(viii) Dynamic time delay for the LST is approximately 15 sec, 

whereas that of the a£ tracker is approximately 6 sec. This 
applies to the transient behavior. 

Figures 51 and 52 show the sample and statistical error time plots 
for the 60 fps climb profile. Initial and final altitudes are 1,600 and 
10,360 ft. It is noted that the simple altitude select/hold logic used 
in the simulation overshoots the reference altitude by 350 ft. This is 
caused by the nonlinear elements (velocity and acceleration authority 
limitors) in the control loop. The following comments and remarks apply 
to this case. 


(xi) By comparing the altitude rate estimates in Fig. 51. the LST 
tracker estimate (in the steady state region) is noisier 

than that of the a$ tracker. This is caused by the deteriora- 
tion in the level switching time computation accuracy. 

(x) The time delay effect (-6 sec) is very close to that of the 

oB tracker except for the "exponential" tailing off when going 
level. The latter is caused by the level occupancy time be- 
coming longer as the altitude rate becomes smaller. 

(xi) The standard deviations of the rate errors were 3.1 and 
4.1 fps for the a$ and LST trackers respectively. The peak 
rate errors were 14.7 and 21.8 fps. The comparable altitude 
errors where 24.8 and 17.4 ft for the a0 tracker and 29.5 
and 67.1 ft for the LST tracker. 

(xii) The transient period for the LST tracker is longer than 
that of the a$ tracker. The average periods were 26.7 and 
32.4 sec. The transient period is defined to be that period 
where there is 6 fps and the larger rate errors. 


Table 11 shows the summary of the statistical errors for the 
four cases 
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Table 11. Statistical Suamary for Different 
Altitude Rate Profiles 


-5 fps 


10 fps 


20 fps 


60 fps 


oS Tracker 

LST Tracker 

22.9 

( 0 ) * 

25.5 

(30.7) 

1.8 

( 3.9) 

1.4 

( 4.9) 

27.0 

( 9.8) 

25.8 

(52.7) 

2.0 

( 4.4) 

1.3 

( 6.3) 

28.9 

(23.5) 

31.2 

(83.9) 

3.0 

(11.9) 

2.4 

(20.8) 

24.8 

(17.4) 

27.5 

(69.1) 

3.1 

(14.7) 

4.1 

(21.8) 


* Altitude standard deviation (ft) Peak mean error (ft) 

Altitude rate standard dev. (fps) Peak mean error (fps) 

















Conclusions 


Based on the simulation results presented, certain conclusions 
can be drawn. The performance of the aB tracker with the fine resolu- 
tion altitude measurements as input is the basis for comparison. 


(a) The low gain ap tracker configuration performs better than 
the high gain configuration. The altitude rate estimate 
is smoother without sacrificing fast response time. 

(b) If the additive noise is correlated (as in this study) 
rather than independent, the basic aB tracker configuration 
needs to be modified. 

(c) The LST tracker performance is very credible considering 
that it must work with the 100 ft Mode C quantization. In 
the low speed ('5 fps) and medium speed (10~20 fps) regimes 
the estimates are smoother than those of the aB tracker. 

The peak mean errors caused by the time delay due to the 100 ft 
quantization are larger. In the high speed regime (-60 fps), 
the LST estimates are somewhat inferior. 

(d) The basic dynamic delay problem with the LST tracker remains 
the fundamental problem. Compared to the aB tracker the 
delay represents an extra *20 sec at 5 fps and -6 sec at 60 fps 
However, this may be a limit which can not be solved by compu- 
tational considerations alone. For example, it may require a 
cross-link of on-board generated altitude rate estimates. 


In the next chapter. Implications of these errors will be discussed 
from the collision avoidance logic point of view. 




V 


COLLISION AVOIDANCE LOGIC 
AND 

CDTI SENSOR IMPLICATIONS 

In the previous three chapters, various configurations and algorithms 
were discussed to generate state estimates of intruder and Own aircraft 
relative kinematics. Performance comparisons were made in terms of 
"raw” estimation error statistics using Monte Carlo simulations. The 
comparative analyses have provided a great deal of useful information as 
to the relative merits. However, the information per se does not provide 
acceptability of a particular estimation algorithm to certain applications. 
The main purpose of this chapter i" to discuss the estimation algorithm per- 
formance in terms of CDTI and TCAS applications. 

The method of linear error analysis is mainly utilized. This method 
is simply stated as follows: If a dependent variable, y, is a function 

of several statistically independent variables, x^, as 


f (x 1 ,x 2 , . . * , x n ) , 


( 80 ) 


then the errors are related within the first order terms as 

6y * ^f i (x 1 ,. . ,x n )6x i , 
i 

where f ± ( - ) is the partial derivative of f with respect to 
standard deviation of y can be computed by the formula 

\h 


( 81 ) 


The 


” f i (x l» x 2 ,,,,x n )o xi 


( 82 ) 
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Sampling Period Scheduling Logic Application 


One of the unique features of the enhanced TCAS II designed by 

Bendix is its ability to tine the interrogation/reply surveillance schedule 

under microprocessor control. This ability is enabled by a relatively 

narrow interrogation transmission beam-width and the beam stabilization 

with respect to Own aircraft orientation. The narrow interrogation 

beam-width, the selective interrogation time scheduling, and several 

levels of whisper/shout power sequencing give it the ability to operate 

2 

in a very high density (up to 0.4 aircraft /nmi ) traffic environment. 

The scheduling logic criterion is given by 

At * min {At , At, , At , 8 sec} , 
s r b z 

where 

At f = time for the range to change 1000 ft. 

At, “ time for the bearing to change 3 deg, and 

D 

At z * time for the altitude to change 230 ft. 

This is not the sole test. The other test is the target threat status. 

If a target is declared a preliminary threat, then the surveillance of that 
target is cycled at 1 sec intervals. 

Range Sampling Scheduling The range scheduling time is defined 
mathematically by 

At r - | 1000/r | , (83) 

where r is the range rate expressed in fps. In re&lity, the schedule time 
is computed based on the estimates, thus 

At r - | 1000/r | . (84) 

The difference between Eqs (83) and (84) is in the precision loss induced 
by the range rate estimate error. If the difference is denoted by 6At , the 
following first order equation is obtained 
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-T 





r 





6At r 


r 



(85) 


where 5r is the estimation error in the range rate estimate. Dividing 
both sides by A t f , Eq (85) says that the positive (negative) percentage 
error in the range sampling schedule time is identical to the negat ive 
(positive) percentage error in range rate estimate. Based on the Range 
Filter Performance Summary Table 7 in Chapter III, the following Table 12 
can be derived. 


From Table 12 the maximum percentage deviation due to the range 
rate estimation error is 22X; thus, it is apparent that the range error 
would not significantly affect the sample scheduling time computation. 

One note of caution should be mentioned here. The range rate errors are 
computed based on the measurements available in one second intervals. 
Therefore, when the measurement cycle is, for example, 3.7 (or 4) sec, then 
the range rate error is expected to be much larger. Even if the errors were 
four times larger, the changes would be less than 1 sec. (The exception 
is that th-»y would be 1.7 sec deviation for the range cB tracker case.) 


Altitude Sampling Schedule Analogous to the range case a similar 
analysis can be performed for the altitude axis. The altitude error 
equation corresponding to the range error Rq (85) is given by 


6t = 
z 


250 



z 


5z 



( 86 ) 


or in terms of percentage 

100 • (6t /At ) X - 
z z 


100 • (6z/z) % 


(87) 


Referring to the altitude tracker statistical performance in Table H, 
it is clear that the altitude sample scheduling time is not affected for 
the low to medium altitude rate (up to 15 fps) cases. For example, a 
+ 100 X increase in altitude rate error represents 30 fps of estimate instead of 
the nominal 15 fps. Accordingly, the computed sample schedule time is 
9.3 sec instead of the true 16.7 sec. In both cases, the sampling period 
would be 8 sec because of the imposed 8 sec minimum. 
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Table 12. Percent Error in Range Sample 
Scheduling Time Computation (sec) 


Encounter 

Nominal 
Range Rate 
& Sampling 
Time 

Non-Aided 

Target- 6 
Own-Aided 

Range 

oB 

Range 

Square 

<*By 

Tail Chase 

160 kt 
3.7 sec 

4.3 (1) 

4.0 

11.4 

5.6 

Route 

Cross- 

ing 

lx 

280 kt 
2.1 sec 

2.4 

1.8 

6.5 

n 

a 

4.4 

3.4 

12.7 

5.0 

Head- 

on 

lx 

340 kt 
1.8 sec 

3.0 

2.0 

5.5 

2.2 

2x 

6.0 

■a 

10.7 

4.7 

4x 

10.9 

8.4 

21.7 

m 

Parallel 

Turn- 

in 

lx 

0 ' 400 kt 
8-1.5 sec 

7.5 (2) 

(18.2) 

3.1 

10.6 

(9) 

7.4 

(14.6) 

2x 

a 

Ha 

5.9 

1 

■99 

10.6 

(14.6) 



(1) Percentage change to Interrogation Scheduling Time. 

(2) Computed based on the relative range rate of 2C0 kt. 
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When the altitude rate is larger than 15 fps, there could be some 
deterioration. Over estimation of altitude rate (ie., the magnitude of the 
estimate is larger than the magnitude of true state) would in fact help 
the estimation process by providing the measurements at a higher than 
nominal sampling rate. Thus, the deterioration occurs only when the 
altitude rate is significantly underestimated. In the 60 fps vertical 
profile case, the nominal sampling period is 4.2 sec. If the rate is 
underestimated by 15 fps as in the tracker case, then the computed 
period is 5.6 sec with the difference of + 1.4 sec. If the rate is 
underestimated by 22.2 fps as in the LST tracker case, the corresponding 
values would be 6.6 sec and + 2.4 sec. Because of the relatively long 

duration of the transient error with that magnitude, the longer- than- 
nominal sampling period would last approximately 15 sec and 20 sec for 
the ag and LST trackers, respectively. (It is not meaningful to discuss 
the ag tracker results in this application except for comparison purposes, 
because Own measurements are available every second.) 


As mentioned previously, the above discussion is somewhat simplistic 
in the sense that the estimation error was generated based on the one 
second sampling interval and not the computed sampling rate. However, 
the results can be considered the expected errors under the best circumstances. 
It is expected that the transient errors of the LST tracker would be sub- 
stantially worse at long sample intervals. 


Collision Avoidance Logic Applications 

The range and altitude are closely monitored by the collision 
avoidance system (CAS) logic in order to ascertain an intruder's threat 
status. Actually the CAS logic consists of two parts. One is to determine 
the presence of a threat (called Threat Detection Logic) , and the other is to 
determine an escape maneuver (called Resolution Determination Logic). The 
main concern here is the detection performance. Basically, the (vertical CAS) 
detection logic consists of two parts. One is the so-called range test, 
and the other is the altitude test. 
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When the range is doling, the range test monitors the range closure 
time called tau (T r ) . Mathematically the test is expressed as 



where d^ is the minimum range guard and 9^ is the threshold value. When 
T r is less than the threshold, 6 t , the intruder is said to pass the 
range test. The value of d f varies from 0.075 to 1.3 nmi, and 0^ varies 
from 18 to 35 sec, both depending on own altitude. 


The main idea of the altitude test is to examine the projected relative 
altitude t sec into the future. If the projection is within a threshold 
then the intruder is said to pass the altitude test. Mathematically, the 
altitude test is expressed as 



Z T^ 


+ T 



(89) 


The altitude threshold (6 ) value varies from 750 to 950 ft, again depending 

£ 

on Own altitude. 


Because these tests utilize the estimates, it is necessary to examine 
the effect of the estimation errors on the tests. These are discussed in 
the following sections using linear error analysis methods. As typical 
parameter values, 1 nmi, 30 sec and 850 ft are used for the minimum range 
guard, the projection threshold and the relative altitude threshold, 
respectively. 


Range Closure Time The linear perturbation error equation corresponding 
to Eq (88) is given by 

6t ■ - — (fir - t fir) . (90) 

r r r 

Assuming that fir and fir are independent fl-.e., E{fir fir) * 0 ) the standard 
deviations are related by 
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{^1 


0 

T 


i [a 2 ♦ , 2 O. 2 ] l ' 2 

• r r r 


(91) 


Referring to Table 7 in Chapter III and using 30 sec for t , the 
following Table 13 is obtained. The error for the parallel turn-in case 
was computed using the root mean square instead of the standard deviation 
values at 200 kt. The following comments and remarks apply. 

(i) All the tracker configurations performed well for rectilinear 
encounter kinematics. The worst was 6.5 sec error for the 

a 6 tracker for the head-on case when the measurement noise 
was four times the nominal. 6.5 sec is 22% of the nominal 
30 sec. Thus, for this case, 34% of the time, the threat 
warning would be delayed 6.5 sec (i.e., there would be 
23.5 sec left for appropriate escape maneuvers.) 

(ii) The closure time errors for the parallel turn-in case vary 
from 4.2 to 6.3 sec for all configurations except the aided 
tracker. Half' of the errors are attributable to the dynamic 
delay error. This implies the errors are sustained for the 
duration of the maneuver. Therefore, the averaging of the 
closure times over a few measurement intervals would not im- 
prove the error. Therefore, this would be a vulnerable period. The 
degradation of the aB tracker performance is less than the other 
two trackers. This is due to the high feedback gains. 


(iii) The closure time errors are computed on a steady state 

basis. The errors would be larger during the initial transient 
period. (This is especially true for the Range square aBy tracker.) 
This implies that the CAS protection is not reliable for 
pop-up targets. Pop-up targets are the ones which remain 
in the antenna s shadow until it is at a very close proximity of 
Own. In this case, the trackers would not have sufficient time 
to settle the initial transient errors. 


In conclusion.lt can be stated that all of the trackers provide 
good range protection for steady state rectilinear encounter kinematics. 
The protection becomes severely deteriorated during maneuvering periods 
except for the aided tracker configuration. Because of the long settling 
time, the Range square aBy tracker does not provide good protection for 
very close pop-ups. 
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Table 13. Range Closure Time Errors (sec) 


1 " " ■ ■ ■ — - 1 

Encounter 

Nominal 
Range Rate 

Non-Aided 

Target- & 
Own-Aided 

Range 

o8 

Range 

Square 

o8y 

Tail Chase 

160 kt 

1.3 

1.2 

3.4 

1.7 

Route 

Cross- 

ing 

lx 

280 kt 

0.7 

0.6 

2.0 

0.8 

2x 

1.4 

1.0 

3.9 

1.5 

Head- 

on 

lx 


0.9 

0.6 

1.6 

0.7 

2x 

340 kt 

1.8 

1.3 

3.2 

1.4 


4x 


3.3 

2.5 

6.5 

2.8 

Parallel 

lx 

0 - 400 
kt 

5.7 

0.9 

4.2 

4.9 

Turn- 

in 

2x 

6.3 

1.8 

6.1 

5.4 


Altitude Projection Error Analysis The altitude protection is 
provided by monitoring the projected relative altitude. The projection 
time is when the range is expected to be the minimum. See Eq (89) . The 
effect of the estimation errors on the altitude protection can be 
analyzed by examining the linearized projection errors. 

Assuming the projection period of 30 sec and utilizing the 
statistical errors for the altitude trackers given by Table 11, the 
following Table 14 of the projection errors is obtained. 

The following remarks and comments are derived from the results. 



Table 14. Altitude Projection Errors (x«30 sec) 


Vertical 

Rate 

aB Tracker 

5 fps 

76.9 (1) (117.0) (2) 

140.0 (3) 

10 fps 

87.0 (141.8) 

166.4 

20 fps 

118.9 (380.5) 

398.6 

60 fps 

117.8 (458.4) 

473.3 


LST Tracker 


67.5 (177.7) 

190.1 


64.8 (241.7) 

250.2 


103.2 (707.9) 

715.4 


152.5 (723.1) 

739.0 


(1) standard deviation 

(2) worst mean error 

(3) root mean square of (1) and (2) 


(i) The range of errors for the aB tracker is 


standard deviation 
peak mean error 
rm9 of these two 


76.9 - 118.9 ft; 

117.0 - 458.4 ft; 

140.0 - 473.3 ft. 


These errors apply to projecting Owr altitude 30 sec into 
the future. 


(ii) The range of errors for the LST tracker is 


standard deviation 
peak mean error 
rms of these two 


64.8 ~ 152.5 ft; 
177.7 - 723.1 ft; 

190.1 - 739.0 ft. 


These errors apply to projecting the intruder altitude 30 sec 
into the future. 


4 

4 









till) If the above sets are combined in an rms sense, then 
-\e relative projection errors are obtained. 

standard deviation: 100.6 ~ 193.4 ft; 
peak mean error : 212.8 - 856.2 ft; 
rms of these two : 235.4 ~ 877.7 ft. 


(iv) The worst combined errors represent 27.7 to 103. 3% of 
the altitude separation threshold of 850 ft. 


(v) The major contributor of the combined errors is the 
peak mean errors of the LSI tracker. The source of 
this error is, of course, the inherent tracker error 
introduced by the 100 ft quantization. Thus, during 
this transient period, the LST estimates would not 
provide sufficient protection. 

(vi) During the constant altitude rate flight (for both 
intruder and Own) , the combined errors are less than 
200 ft (23.5% of 850 ft). If no other error source 

is present, the estimation precision may be sufficient 
to gauge the threat situation. 


(vii) The estimation error statistical data were obtained by 
assuming a high frequency error of + 23.8 ft + la 
and the Mode C 100 ft quantization. Low frequency errors 
such as bias, scale factor or pressure transducer 
dynamic delay errors were not considered. For complete 
assessment, these error sources need to be factored into 
the analysis. 


Three conclusions can be drawn based on the above simple analysis. 

(a) If both intruder and Own maintain steady altitude rates, 
then the vertical threat assessment can be made with 
sufficient precision. 

(b) During transient periods, the dynamic delay errors may not 
allow accurate threat assessment . 

(c) If the combined low frequency (bias, scale factor, or drift) 
error is 200 ft or larger, then the vertical threat assessment 
accuracy becomes marginal. 
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One vs. Two Dimensional Horizontal Test The range test given by 
Eq (88) is called the modified tau test . The more standard test is 
given by 


T 

r 




(92) 


x f can be Interpreted to be the time to collision (r + x^ r*0) . When 
the relative kinematics are non-accelerating and when the aircraft are in 
collision courses, the above interpretation is correct. However, it 
is no longer correct if the miss-distance is non-zero. In fact, the test 
Eq (92) is inefficient in the sense that it passes many non-threatening 
intruders; these are unnecessary alarms. 


The one dimensional range test Eq (92) is designed for range-altitude 
TCAS, especially the so-called active BCAS which is the direct predecessor 
to the current minimum TCAS II. It does not require a directional (bearing) 
capability. In light of the added enhanced TCAS II capabilities, a 
better two dimensional test is available utilizing accurate bearing 
measurements as well as stabilization with respect to Own attitude 
orientation. 


Now assuming rectilinear motions, range is given by 

r(t) - l(x + tx) 2 + (y + ty) 2 ] ^ 2 , (93) 


where x, y, x and y (A notation is dropped here) are constant. The 
so-called time to closest point of approach 1® defined to be the 

value of t which minimize r(t). The miss distance (m^) is the minimum 
range. These two quantities are computed by the following formula: 

t cpa “ + yy > /y2 ; m d " l*y - yx|/ v ( 94 > 


where v is the relative speed defined by 


- [x‘ + 


* 2 , 

y J 


1/2 


It is noted that x^p^ and x f are in general different. 
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The two dimensional horizontal test would take the following form 


T_ , < 0 

CPA — t 


.AHD. 


(95) 


where 0 T is the same as in Eq (88) and (92), and 6 m is the minimum 
miss-distance threshold. (0 m may or may not be the same as d f .) 
Actually, the condition "m d <_ 0^" is sufficient to identify threat 
status as long as is positive. The condition " T cp^ 1 ® T " states 

that one can wait until © T sec to go before he must take an evasive 
action. 


The range can be expressed in terms of t cpA and as 
r(t) - [m* + v 2 (t - t^) 2 ] 1/2 , 


(96) 


by defining time-to-go, t G - - t , then range-to-go becomes 


z . r 2 . 2 .2, 1/2 

r(t G ) " l “d + V *G 


(97) 


The range closure time, T r , is given by 


2 2 2 
m d + v C G 

v 2 t„ 


(m d /v) 2 + t 2 


(98) 


Figure 53 shows the relationship between t r and t^. When the miss distance 
is zero (i.e., collision), then the range closure time is the same as 
time-to-go. 


When an intruder is identified to be a near-miss by the two dimen- 
sional test Eq (95), then it automatically satisfies the modified range tau 
test Eq (90), if the intruder is closing at all. However, the standard 
range tau test (92) may or may not be satisfied. In the case of a near-miss 
with m d ■ 0.3 nmi and 0 £ - 30 sec, the test is satisfied if the relative 
speed, v, is equal to or greater than 72 kt; otherwise, it would not be. 
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m d /v 


Time-to-go 


Figure 53. Range Closure Time as a Function of Time-to-go 

There are numerous (infinite) unnecessary alarm examples of passing 
the standard or modified range test but not the two dimensional test. A 
typical example is shown in Fig. 54. At t-0, Own flies due north at 
200 kt. Intruder flies with 150 deg heading at 300 kt. The initial intruder 
position is exactly 3.5 nmi due north of Own. The following CAS para- 


meters can be computed: 



miss distance, m^ 

as 

1.08 nmi 

time to CPA, 

- 

24.75 sec 

range @ t-0 

■ 

3.5 nmi 

range rate @ t-0 

- 

-459.8 kt 

standard tau @ t-0 

m 

27.38 sec 

modified tau @ t-0 

m 

19.56 sec. 


If the miss distance threshold is 1 nmi, the intruder will not be 
classified as a threat, since « d > 1 nmi. t CpA states the miss distance 
is reached in 25 sec. However, if the tau threshold is 30 sec, then both 
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l*-V 


x (North) 
t 
I 



3.5 nmi 


200 kt 


t 


y (East) 


Own 


Figure 54. Relative CAS Geometry for the Counter Example 
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the standard and modified range tau tests are passed, i.e., this intruder 
will be classified as a threat. Therefore, this is an instance of 
unnecessary alarm. The above discussion and example show that the two 
dimensional horizontal test is a stronger and a better threat detection 
criterion than the one dimensional range tau tests. 

Error Analysis The next step is to determine how well these para- 
meters can be determined from the state estimates. The linear error 
analyses of these parameters with respect to the estimation errors becomes 
very complex. Therefore, the Monte Carlo simulation was performed. The 
simulation scenario wa9 the head-on encounter (Case (C) of Fig. 14) . True 
miss-distance is 0.31 nmi. 

The Monte Carlo simulation program was modified to compute the miss- 
distance (m d ) and the time to the Closest Point of Approach (t cpa >- These 
are computed based on the true state variables, the non— aided x— y filter 
state estimates and the Target— and Own-Aided x— y Ka lm a n filter state 
estimates. The measurement update period is 1 sec. Simulations were run 
for the nominal, twice the nominal and four times the nominal range and 
bearing errors. To repeat, the nominal is defined to have a ranging error of 
+ 75 ft (+lo) and a bearing error of + 1 deg (+lo) . The errors are 
assumed to be Independent white noise. 


Table 15 shows the summary of simulation results. It is organized 
to show the error dependence on the range as well as the noise level. For 
example, at range of 5.7 nmi, the time to the closest point of approach is 
60 sec. At this point in this particular encounter, the m d and t cpA errors 
based on the Non-Aided filter estimates with the nominal measurement errors 
were 0.35 nmi and 1.87 sec, respectively. This means that 68Z of the time 
the m d estimate is between 0.0 to 0.66 nmi compared to the true value of 
the 0.31 nmi. If the miss-distance protection is 1 nmi, then 951 (2o band) 
of the time the test would make a correct assessment. The following 
consents and remarks are derived from the simulation results. 
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V 


Range 

Non-Aide< 

Filter 

Aided 

Filter 

Measurement 

t cpa 

“d 

(nmi) 

t cpa 

(sec) 

“d 

(nmi) 

!cpa 

(sec) 

Error Level 

5.7 nmi 

0.35 (1) 

1.87 (2) 

0.24 

1.13 

nominal 

0.73 

4.17 

0.43 

2.40 

2 times 

60 sec 

1.48 

9.53 

0.90 

5.33 

4 times 

4.2 nmi 

0.27 

1.13 

0.19 

0.67 

nominal 

0.51 

2.67 

0.31 

1.43 

2 times 

45 sec 

1.03 

6.6 

0.64 

3.53 

4 times 

2 . 7 nmi 
30 sec 

0.16 

0.73 



0.11 

0.43 

nominal 

0.28 

1.53 

0.20 

0.90 

2 times 

0.52 

3.47 

0.34 

2.03 

4 times 




















































(i) The errors are, for all practical purposes, proportional 
to the TCAS sensor error magnitude. 

(ii) At the range of 4.2 nmi (45 sec to CPA), the non-aided 
tracker would make a 100% correct assessment , if the 
measurement errors are nominal. This drops to 68% at twice 
the nominal. At four times the nominal values, the ratio 
would be less than 50%. 

At the same range, the aided filter would make a correct 
threat assessment of 100% up to two times the nominal. 

At four times level, the reliability drops to 68%. 

The T rt> . errors are not significant, the maximum being 
CPA 

6.6 sec. This would still afford sufficient protection 
time. 

(iv) At the range of 2.7 nmi (30 sec to CPA), the estimates 

improve substantially. Both trackers would make a correct 
assessment of almost 100% of the time with up to two 
times the nominal error. 

At four times the nominal, the reliability drops to '.*5% 
and 90% for the non-aided and aided filter respectively. 

At this range, the worst error was 3.5 sec for the non- 

aided tracker at four times the nominal measurement error 
level. The error represents 11%. This magnitude of time 
error is thought to be noncritical. 


It should be strongly emphasized that the above comments are based 
on a single simulation. The reliability numbers are based on the assump- 
tion that the errors arc normally distributed. The analysis should not 
be taken at the 10*^ type precision. With these caveats, an important 
conclusion emerges; at the critical time of 30 sec-to-go, the two dimen- 
sional tests can assess threat status fairly accurately. One corollary 
is that the threat assessment becomes more accurate as the range becomes 
closer. 
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CDTI Applications 


The CDTI's major function is to provide the surrounding traffic informa- 
tion base to the pilots. Therefore, the CDTI estimation accuracy require- 
ments depend on how the pilots choose to use it. Obviously very accurate 
position information is not needed if the pilots want to be simply aware of 
proximite traffic. If the CDTI is used to perform sophisticated Electronic 
Flight Rule (EFR) tasks, then highly accurate estimates may be needed. The 
EFR tasks may include self-spacing along a route, route-crossing, merge into 
a traffic stream and so on. 

Because of the human factor element in the CDTI applications the only 
meaningful way to determine the CDTI accuracy requirements is through pilot- 
in-the-loop simulation study. However, in this environment, other issues 
come into play such as display size, brightness, contrast, symbology and 
so forth. A recent simulation study [ 20 ] obtained two results concerning 
parameters relating to CDTI estimation accuracy requirements. 

(a) Displayed traffic position errors with standard-deviation values 
up to 0.3 nmi range and 8° azimuth had negligible effect on the 
ability of the pilots to perform the self-spacing task; and 

(b) Display of the lead aircraft groundspeed was found to affect 
the mean spacing performance, especially during periods of 
speed or spacing changes. Pilot comments cited the ground- 
speed information as a definite aid in performing the spacing 
task. 

These must be considered as two separate factors, since it is not possible 
to obtain accurate ground speed estimates when the measurement error 
magnitudes are 0.2 nmi and 8 deg for the range and bearing, respectively. 

It is also interesting to note that the ground speed information was used 
by the experiment pilots as a damping signal to prevent spacing "overshoot". 

It is safe to say that the accuracy per se will not be a major issue 
unless, for example, an in-trail following "flight director" signal was 
generated and incorporated as an integral part of the CDTI symbology. In such 
a case, the connections between the (pilot) performance, the flight 
director accuracy, and the underlying state estimation accuracy can be 
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analyzed in a more qualitative fashion rather than through a more 
quantitative simulation analysis. It 1s very difficult to specify 
estimation accuracy requirements for the CDTI applications in a com- 
prehensive quantitative way. 

In the following sections, effects of measurement and corresponding 
estimation errors are discussed with respect to selected CDTI variables. 

Also discussed are smoothing (as opposed to filtering) algorithms for 
selected CDTI applications. 

Bearing Errors In many cases, the surrounding traffic Information 
is superimposed with other symbology on an EHSI (Electronic Horizontal 
Situation Indicator) or an MFC (Multifunction Display). Other information 
may include reference air-route, map and terrain information, waypoints 
and navaids. Thus, the traffic position would be referenced to 
a local- level map fixed coordinate system. (See Fig. 55 .) This can be 
accomplished by transforming the TCAS measurements - relative range (r) , 
bearing (b) and altitude (z) - to a north referenced Own fixed local level 
coordinate system utilizing Own body attitude angles. To this relative 
NED position is added Own earth-fixed NED position to obtain the "true" 
horizontal projection. This process is explained in Appendix A. (See Fig. 
A-l.) When Own attitude angles are not properly accounted for, then the 
horizontal projection could be substantially in error [21 j . The angle, b-HJ/ 
(TCAS relative bearing plus Own heading) , is not the true horizontal 
bearing with respect to north, when the Own roll or pitch angles are 
non-zero. For a 20 deg roll angle, for example, the peak error (depends 
on (<Jj) would be 5 deg if target elevation is 10 deg (Fig. 56) • Thus, it 
is safe to say that proper transformation must be performed either within 
or outside the TCAS processor for the CDTI application. 

Assuming that the above problem Is solved, we need to discuss high 
frequency error magnitude. The draft TCAS MOPS specifies the error magnitude 
as 9 deg rms. Flight test results of one pre-production TCAS, designed 
by Dalmo Victor, showed the bearing error of 5-10 deg rms [22]. Bench test 
results of TCAS Engineering Unit designed by Bendix showed the error 
magnitude between 0.6-2 deg (lo) [23]. Some of these numbers are preliminary. 
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Considering the findings of the aforementioned NASA study, these error 
magnitudes are within the maximum allowable for the CDTI applications, 
provided that some filtering is performed prior to displaying it to the 
pilots. Certainly, these are sufficiently accurate in terms of clock angle 
indications to aid pilots in VFR applications. 


According to the TCAS MOPS, the following Cartesian ct8 estimation 
algorithm has been used successfully to develop smoothed target bearing 
estimates from reply measurements. 



i 


Initiation The bearing and range measurement of each reply are 
used to form x and y position: x = r cos(b) and y * r sin(b) . 

The first three measurements are used to form a least-squares 
estimate of the x and y positions and velocities. 

Prediction A predicted x,y position for the next scan is formed 
by adding the product of the last-scan velocity estimate and the 
time since the last scan to the last-scan position estimate. 

Update If a valid measurement: is available it is combined with 
the range measurement to form the x and y measurements. The update 
is made using a standard a& tracker algorithm in both x and y. 

The gains are the same in the x and y coordinates. The gains vary 
during the first eight scans following initiation, and then are held 


constant, as shown below. 



Track age 

a 

6 

4 sec 

0.700 

0.300 

5 

0.600 

0.200 

6 

0.524 

0.143 

7 

0.464 

0.107 

8 

0.416 

0.083 

9 

0.378 

0.067 

10 

0.345 

0.054 

11 or more 

0.318 

0.045 

The position update equation 

in x is: 



x(t) estimate * x(t) prediction + [o * (x(t) measurement - 
x(t) prediction)] . 

The velocity update equation in x is: 

x(t) estimate - x(t-T) estimate + l -| * (x(t) measurement - 
x(t) prediction)] , 
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where T is the time difference between the current and previous 
measurements. 

The y equations are analogous. 

Bearing reputation The updated x and y are converted to bearing 
using the tangent function. 

It should be commented that the above filtering algorithm is 
used to correlate transponder replies and target tracks Internally stored 
in the TCAS surveillance processor. The algorithm would not provide good 
position (or bearing) estimates for the CDTI applications if Own aircraft 
undergoes a maneuver (either pitch or roll) . This point was explained 
previously in this section. A better way to obtain required estimates for 
the CDTI applications is (1) to obtain the north referenced local level x 
and y measurements using Own. <|>,6and i|> and r and b from the TCAS , then (2) 
to use the aB tracker algorithm described in Chapter II to obtain the hori- 
zontal state estimates from which the bearing estimate can be computed. As 
mentioned previously, the attitude stabilization can be performed in the CDTI 
processor module outside the TCAS as long as timing is synchronized. 

Target Prediction Vector Error The target prediction vector is 
sometimes included in the CDTI symbology - see Fig. 55. The Own pilot 
can extract valuable target short term future information from this vector. 

The prediction vector is computed based on target ground speed (vq) 
and ground track angle (i|> G ) . These variables are given by 

0 G - Kx 0 + Ax ) 2 + (y Q + Ay) 2 ] ^ 2 , and (99) 

<J> G - tan -1 [(y Q + Ay)/(x Q + Ax) ] , (100) 

s s 

where x n and y. are Own horizontal velocity components provided by the on- 
board navigation system, and Ax and Ay are the TCAS relative velocity 
estimates. Obviously* the ground speed and ground track errors would depend 
not only on the TCAS errors but also on the on-board navigation system. 
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The navigation errors vary according to available navaids, configuration 
and geometry, body rate sensor or INS, and so on. The target ground 
speed is sometimes included in the target data tag. Numerical values of 
the ground speed and flight level are shown along with the target 
identification. Thus, the Own pilot would have a ready reference with 
respect to his Own instrument reading. 

In order to assess the TCAS sensor effect on the above parameters, a 
simple in-trail following scenario was incorporated into the Monte Carlo 
simulation program. The ground speed and track angle estimates were 
computed according to Eqs (99) and (100) , except the Own velocity 
components were used. The error statistics between the estimated and true 
variables were obtained from this simulation. 

The in-trail following scenario was a simple one simulating downwind 
turn to final. The lead aircraft was placed 3.8 nmi ahead of Own flying 
due north at 200 kt. After flying straight for 1.8 nmi, the lead executes 
a 170°, 15° bank angle left turn, and then flies due south. Own follows the 
lead by flying due north at 220 kt and executing a similar turn at approxi- 
mately the same location. These were performed open-loop, and the initial 
conditions and flight parameters were made slightly different so that 
the lead and Own would not traverse the same trajectory. Figure 57 shows 
the horizontal projections. 

Figure 58 shows statistical error time plots of the ground speed 
and track angle estimates based on the non-aided x y tracker configura- 
tion. Only the non-aided tracker results are presented and discussed. 

For the aided tracker configuration, the estimation accuracy of the tar- 
get ground speed and tracker angle does not directly reflect the CDT1 
position sens. r accuracy. 
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Figure 57. CDTI In-Trail Following Scenario 



IP 





The simulation data were taken for three levels of sensor measurement 
grfor* nominal, twice, and four times the nominal* The following 
comments and remarks are derived from the results. 


(i) As expected, the high frequency errors are dependent on 
the relative position and velocity. When the lead is 
directly ahead of Own and aligned, the ground-speed 
error is proportional to range error, but the heading 
error is more affected by the bearing error. However, 
when the alignment becomes more oblique, the range and 
bearing errors are no longer orthogonal; thus, mixing 
into both speed and angle estimates. 

(ii) This can be seen in all three ground speed error plots. 

After the initial transient period, the range standard 
deviation is sr *11; as the bearing deviates from 0 as the 
lead turns, the error becomes larger; finally, as both 
aircraft line up in the same direction after the turn 
the error becomes smaller again. The opposite behavior 
is true for the heading error. After the initial 
transient period, the heading error is relatively large; 
as the bearing becomes more oblique, the heading error 
becomes smaller; and finally, when both aircraft line up, 

it becomes larger. 

(ii) The standard deviations for the ground speed were 13.0, 

21.5, and 41.7 knots, and the average peak mean errors 
were 19.4, 21.3, and 30 knots for nominal, twice and 
four times error levels, respectively. 

The standard deviations for the heading were 5.8, 10.3, 
and 19.9°, and the peak mean errors were 5.4, 6.3, and 
8.7°. 

The standard deviations of the heading error are roughly 

proportional to the sensor-error level. The ground 

speed error shows a trend, but it is not as clear-cut 

as for the heading errors. The peak errors are similar; for the 

nominal and twice nominal cases. They are 19.4 vs 21.3 knots 

and 5.4 vs 6.3°, but at four times the nominal, the peak 

mean errors increases almost 50%. 

(iii) If the numerical value of the target ground speed is to 
be shown as a part of target-data tag, then it may be 
advisable to quantize the magnitude so that pilots are 
not annoyed by the noisy digital indication. 
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(iv) The third row of the summary Table 16 shows the 60 sec 

projected position as based on 200 knot target speed. The 
table shows the standard deviation and the maximum rms magni- 
tude. Accordingly, for example, in the case of the nominal 
sensor error level, the tip of the target prediction vector 
would, on the average, move +0.4 nmi (+lo) from one sample 
period to another. This implies that, at a display scale of 
1 nmi/inch, the vector tip hops around 0.4 inch. 

This indicates that the digital display of the ground 
speed and the target prediction vector may be pilot 
selectable CDT1 functions, i.e., the pilot might choose 
to suppress the display elements if the signal quality 
becomes below his acceptable level. 


Table 16. Summary of Prediction Vector Error 
for Non-Aided x y Tracker 



Nominal 

Twice Nominal 

4 Times Nominal 

Ground-speed Error 
knots 

13.0^^ 
(19! 4) (2 ^ 

21.5 

(21.3) 

41.7 

(30.0) 

Heading Error 
degrees 

5.8 

1 

10.3 

19.9 

Prediction Position 
Error (60 sec) nmi 

0.4 1 2 (3) 4 

<0.49) <4) 

0.70 

(0.86) 

1.35 

(1.53) 


(1) Average standard deviation. 

(2) Average peak mean error. 

(3) Standard deviation of 60 sec prediction position error. 

(4) Maximum rms 60 sec prediction position error. 

The following are tentative conclusions. It should be noted that 
they are based on simulation results of one particular scenario, and the 
TCAS/CDTI sensor is the only error source. 

(a) The target ground speed and heading angle estimates may 
provide useful CDTI information if the sensor error 
magnitudes are less than + 150 ft (+lo) and +2* (+ic) for 
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the range and bearing, and if the underlying relative 
geometry is favorable. 

The target ground speed and heading-related CDTI display 
parameters, such as the ground-speed data tag or the 
target prediction vector, may (should) be made pilot 
selectable functions. Pilots would be able to decide 
at what noise level these signals cease to become use- 
ful for their tasks. 

The conslusions are supported to a certain extent by a simulation study. 

A NASA study ( 24] showed that the pilot's intrail-following performance 
did not deteriorate much (qualitatively speaking) for up to 20 knots of 
target ground-speed indication. 

Application of Smoothing Algorithms . In some CDTI applications, 
the state estimates at some, past time are more important than those at 
the current time. This is the case with the so-called Constant Time 
Delay in-trail following task. Very briefly, this criterion states that 
the Own follows the Lead T seconds later. Mathematically, it can be 

expressed as: • 

I 

d Q (t) = d T (t-T D ) (101) 1 

where dg(.) and d^.(.) are the distances traveled by Own and the Lead 
along a fixed air route, and T^ is a fixed delay time. The current Own 

position is where target was T^ sec ago. One of the major advantages ^ 

of this criterion compared to others (for example, a Constant Time Predictor) 
is that the velocity profile must be identical except for the time delay, 
i.e. , : 



v Q (t) “ v T (t-T D ) • 


( 102 ) 
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Reference ( 25 ] discusses analytical aspects of the various in-trail 
following criterion. 

From the regulator control viewpoint, the perturbation acceleration 
command, 6a c , to satisfy the criterion of Eq (101) may be expressed as: 

5 a c<t) * K p [d 0 (t)-d T (t-T D )]+K D [v 0 (t)-v T (t-T D )] (1C 

where Kp and are the proportional and derivative regulator gains. 

Equation (103) may be a basis for a flight director design via a 
speed-error tape indicator, or an outer-loop guidance law design; how- 
ever, it is not advocated that Eq (103) is implemented exactly. Exact 
implementation would depend on other factors such as the inner-loop 
design . 

It is interesting to obtain a rough idea of sensor noise effect on 
the acceleration command error. The following are assumed: 

position error = 60 ft (lo); 
speed error = 13 knots (lo); 

2 “2 

proportional gain (K ) = (0.2) sec ; and 

P -1 

derivative gain (K^) = 0.2 sec 

Then the acceleration command error standard deviation, a is 

CL 

given by: 

o a - [Kp 2 0 d 2 + - 2.6 knots/sec (1 

Obviously, this would be excessive. Most of the above error is 
attributable to the velocity error. 
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In Che majority of the CDTI symbology devised for simulation 
studies, the pilots extract the necessary control parameters from the 
history dots with one exception [26]. In the exception case, a flight 
director type display was also used in conjunction with the history 
dots. The Lead past positions are displayed as d^ts. The one 
corresponding to T^sec ago is marked specially, and the target ground 
speed alphanumerics may be shown by the special symbol as the speed 
reference. The pilot tries to put Own aircraft symbol on the dot with 
the proper speed to satisfy the criterion. 

Those display parameters (especially the ground speed) were generated 
based on simulated true values. However, in actuality, these must be based 
on the estimates. There are essentially two methods of generating the position 
and velocity estimates, d T (t-T D > and v^,(t-T D > . One is to use the 
filtering algorithm, i.e., 

E(d T (t-T D ) | measurements up to time t-Tp). (105) 

The other is to use the smoothing algorithm, i.e., 

E{d T (t-T D j | measurements up to time t} (106) 

Here, the notiation E{ } means the standard conditional expectation. A smoothing 
algorithm of particular interest is the so-called fixed time lag, fixed interval 
smoother [27 ]. The fixed time lag means that only the estimates at time 
t-Tp are computed as the current time, t, advances. The fixed interval means 
that the data interval is fixed rather than extending all the way to the 
initial time. Usually the fixed interval is taken symmetric with respect to 
the reference time, t-Tp. 


154 


A 


-nr 


ORIGINAL PAGE .5 
QE POOR QUALITY 




Appendix C derives a proposed recursive fixed time lag, fixed inter 
val smoothing algorithm based on a linear three states Newtonian dynamic 
model. Appendix D derives a proposed nonlinear smoother algorithm based 
on a circular-arc trajectory dynamics model. It incoporates a parameter 
identification subalgorithm based on the bank-of-Kalman-f liters idea. 
Readers are referred to the Appendices for detail. 


During the course of this study, these algorithms were implemented 
to test their applicability. The following comments pertain to this effort. 


(i) The recursive algorithm experienced a numerical stability 
problem. One cause of this problem was thought to 

be the fact that the closed form poles (i.e., eigen- 
values of the system matrix, F, in Eq. (C.27) are all 1. 
However, a semirecursive least-square smoother algorithm 
based on Eqs (C.13), (C.23), and (C.8) did not experience 
numerical instability problems. In this approach, the 3x3 
matrix inverse may be precomputed and stored. 

(ii) Because of the modeling errors during the turn maneuver, 
the usable smoothing interval was 7-11 sec, (i.e., half 
intervals of 3-5 sec at the nominal error levels). With 
these short smoothing intervals, the estimates were not 
much more accurate than the fixed gain filtering algorithm. 

(iii) Two methods to use nonlinear interpolation algorithms were 
attempted as described in Appendix D. One was based on 
Fig. D-l, and the other was a linear least squares id- 
entification algorithm based on Eq (D.ll). 

(iv) Because of the dynamic delay problem associated with the 
step change in turn-rate, to, the smoothing interval was 
short compared to the linear case. 

(v) Major problems associated with the bank-of-Kalman-f liters 
decision process (Fig. D-l) were that the individual rms 
error distribution had more than one local minima. Further 
more, these minimum rms errors showed substantial random 
nature from sample to sample. 

(vi) The effort to identify the unknown coefficient (cos wA) in 
the auto regressive equation (D.ll) using a standard 
technique (28] was not effective. The sensor did not 
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provide sufficient precision to identify the value of 
cos o>A » 0*9986 for u> * 3°/sec. In order to obtain the 
statistical precision of this order, one would have to 
use thousands of measurements during which a> would have 
changed . 

The basic conclusions from the smoother study are two: (1) smoothing 

algorithms are of limited value for the CDTI applications, and (2) it is not 
possible to estimate (filtering or smoothing) the turn rate from the TCAS 
measurements. Compared to the previous results on the position and velocity 
estimates based on the Mode S ground sensor [3] , the possibility of obtaining 
usable turn rate estimate is very small for the TCAS sensor. The differences 
are 0.04 vs 1° bearing error magnitudes and 4.6 vs 1 sec sampling times for 
the ground-based system and the TCAS. 

From the data storage point of view, to generate the past T^ = 60 sec 
history dots at 4 sec apart, the smoother algorithm required substantially 
more memory — at least 280 cells compared to 60 cells for the filtering algorithm. 
For the filtering algorithm, the target x y position and velocity estimates (with 
respect to the underlying map) need to be stored every 4 sec. For the smoothing 
algorithm (assuming a 10 dots ahead and 10 behind algorithm) , the last 70 
sec worth of relative range and bearing measurements as well as pseudo x y 
measurements need to be stored. The relative range and bearing are used to 
compute measurement error covariances* The computational requirement for 
implementing the smoother algorithm is an additional load. 

Concluding Remarks 

In general, the CAS logic requires a higher accuracy in state estimates 
as compared to the CDTI applications. The difference is functional. TCAS 
must generate an advisory the pilot may follow based on the available 
information, whereas CDTI provides an information base to the pilot 
so that he can make decisions, thus, the former is more tactical whereas 
the latter is more strategic* As such, the TCAS function can be gauged in 
a quantitative manner— for example, the threat detection reliability can be 
related to the estimation accuracy. On the other hand, the measure of merits 
for the CDTI functions are from the pilot utility viewpoint. Thus, for the 
CDTI applications, the estimation accuracy analysis is relative (in the sense 




of human factors research) rather than absolute. It is the subjective opinion 
of this author that a reasonable and accurate traffic sensor (say, one with 
less than 6° bearing error) would provide an accurate enough strategic 
information base to the pilot. 
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VI 

CONCLUSIONS 


In this report various TCAS sensor estimation problems were examined 
from the viewpoint of CDTI and CAS applications. The enhanced TCAS II de- 
signed by Bendix was used as the traffic sensor basis; however, parts of the 
analysis are also applicable to the minimum TCAS II traffic sensor. Three 
problem areas were investigated - horizontal x-y, range and altitude estimation. 
The insight obtained from this study is summarized below. 

In Chapter II, horizontal x and y estimation algorithms were developed 
and analyzed concerning two main factors - signal configuration and filter gain 
selection. By signal configuration is meant any additional information available 
for complementing the basic TCAS measurement of range and bearing. Own and target 
accelerations (or differential velocities) were selected as generic complementa- 
tion signals. The target acceleration signals were assumed to be cross-linked 
via Mode S data link. Three filter configurations were developed based on three 
combinations of different signals. These were (1) non-aided; (2) Own data- 
aided; and (3) Own and Target data-aided. 

Three gain selection methods were developed and discussed for each con- 
figuration - fixed gains, Kalman filter gains, and table-look-up gains* (The 
current Bendix TCAS uses the fixed gain, non-aided configuration algorithm.) 
Performance analysis data were obtained with respect to TCAS sensor noise level, 
TCAS surveillance interval, and Own and Target maneuvers using Monte Carlo 
simulation method. Based on the "raw" error performance statistics, the following 
conclusions are appropriate: 

(i) Combination of the Own and Target data-aided configuration and 
the Kalman gain updating exhibited the best results; 

(ii) Own data-aiding helped when the relative accelerations were due 
to Own maneuver; however, this was not always true when Target 
and Own maneuvered ; 

(iii) Non-aided configuration had good performance when the underlying 
kinematics were rectilinear. It developed large and sustained 
velocity errors, if Own or Target maneuvered; and 
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(iv) Monitoring the measurement residuals for sustained (bias) 
errors is necessary to know when not to use the estimates. 

Also, see Conclusion (vi) below. 

In Chapter III, several range filter algorithms were analyzed. These 
are range and range rate estimated based on (1) non-aided horizontal filter; 
(2) Own and Target data-aided horizontal filter; (3) two state aB tracker; 
and (4) three state range square aBy tracker. Performance analysis data 
were obtained for each algorithm with respect to TCAS sensor noise levels 
and encounter geometries using the Monte Carlo simulation method. Four 
scenarios - were selected to simulate collision encounters: tail chase, route 

crossing, and head-on and parallel turn-in. Based on the "raw" est ima tion 
performance data, the following conclusions are appropriate: 

(v) The aided configuration showed the best performance in all 
cases. The estimation errors were similar in magnitude re- 
gardless of the encounter geometry. The errors were pro- 
portional to the measurement error levels: 

(vi) The non-aided configurations were better than the other two 
range only estimators (the range aB and the range square apy 
trackers) except the parallel turn-in encounter. The diffi- 
culty was the large and sustained mean range rate error 
caused by the filter dyramic delay; 

(vii) The range square filter performed credibly for rectilinear 
encounters. However, this may be due to the gains being 
too low; and 

(viii) The range aB filter suffered from the nonlinear behaviour 
of range at or near the minimum range even for rectilinear 
encounter cases. It exhibited very quick recoveries due 
to its high gain nature; however, due to the very same nature, 
it passed a large portion of the high frequency noise 
achieving less signal smoothing than other filters. 

In Chapter IV, the vertical estimation problems were addressed. The 
main concern was the treatment of the 100 ft quantized Mode C altitude 
reports. A new algorithm called the Level Switching Time (LST) filter was 
designed for this purpose. It was investigated extensively comparing its 
results with those of an aB tracker with non-quant ized altitude input using 
Monte Carlo simulation. The aB tracker is used to obtain Own altitude and 
altitude rate estima < . Thus, it represents the best performance possible 
without complementing it with other signals duch as vertical acceleration. 

The following conclusions apply to this chapter: 
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(ix) If the altitude additive noise is correlated rather than 
independent, the basic a 6 tracker algorithm needs to be 
modified; 

(x) The LST tracker performance is very credible considering that it must 
work with the 100 ft Mode C quantization. In a low to medium rate 

(5 - 20 fps) regime the steady state LST estimates are smoother. 

In the high speed regime (* 60 fps) , the LST estimates are some- 
what inferior. The peak mean errors caused by the time delay due 
to the 100 ft quantization are larger; 

(xi) The basic dynamic delay problem with the LST tracker remains 
the fundamental problem. Compared to the oB tracker, the delay 
represents an extra * 20 sec at 5 fps and ~ 6 sec at 60 fps. 

However, this may be a limit which cannot be solved by algorithmic, 
considerations alone. For example, it may require a cross-link or 
other signals such as the vertical acceleration. 


Various estimates are used to calculate dynamic parameters for other ap- 
plications. In Chapter V, the "raw” estimation performance statistics were 
used to infer impacts to selected CAS and CDTI applications, 
include: 


These applications 


f 

r 


Surveillance Function 

range sampling schedule; 
altitude sampling schedule; 

Collision Avoidance Logic 

range closure time; 
relative altitude projection; 

one vs two dimensional horizontal threat detection; 

Cockpit Display Applications 

bearing errors; 

target prediction vector. 

Both linear error analysis and Monte Carlo simulation methods are used to 
draw technical conclusions. Additionally, smoothing (as opposed to filtering) 
algorithms were investigated for an active CDTI mode application. The algorithms 
are developed in Appendices C and D. The following conclusions are tentative 
to the extent that the statistical base is limited. 
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(xii) 


The next surveillance schedule time computations based on 
either range rate or altitude rate (LST) estimates were 
affected very little. The maximum deviation was 22% for 
the range a£5 tracker results with the four times the nominal 
noise level. However, the results are applicable to esti- 
mates generated with 1 sec sampling interval only; 

(xiii) For the range closure time (tau) computation, all the 
tracker configurations performed well for rectilinear 
encounter kinematics. The worst was 6.5 sec (22% of 30 
sec protection time) error for the ctB tracker for the 
head-on encounter case when the measurement noise was 
four times the nominal. For this case, 34% of the time, 
threat warning would be delayed 6 . 5 sec . 

The closure time errors were 4.2 to 6.3 sec for the parallel 
turn-in case for all except the aided tracker configuration. 

"Half" of the errors were attributable to the dynamic delay. 
Therefore, this error would be sustained for the duration of 
the maneuver; 

The closure time errors would be larger during the initial 
transient period. This implies that the CAS protection 
would not be reliable for pop-ups; 

(xiv) The error in projected relative altitude 30 sec into future 
could be less than 200 ft for steady climb rates for Own and 
target and between 240 to 880 ft during the maneuver transients. 
Compared to the altitude separation threshold of 850 ft, the 
former error magnitude would be satisfactory and the latter 
would not; 

If the combined (Own and target) low frequency error (bias, 
scale factor, or drift) is 200 ft or larger, then the vertical 
threat assessment accuracy becomes marginal; 

(xv) Two dimensional (x and y) threat assessment test was found to 
be superior compared to the one dimensional (range only) test; 
for the head-on encounter case, the threat assessment using the 
horizontal miss-distance was 100% accurate at the range of 2.7 
nmi (30 sec to CPA) for the nominal measurement noise level. The 
reliability dropped to 65% for the non-aided tracker if the 
noise level was four times the nominal; 

The threat assessment became more accurate as the range became 
closer; 
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(xvi) For the CDTI applications, relative bearing should be computed 
within the stabilized local level reference frame. If this is 
not done with the TGAS sensor, this should be done using an 
onboard (navigation) computer. The error could be as large 
as 5 deg for the Own roll attitude of 20 deg; 

(xvii) The non-aided configuration velocity errors in terms of ground 
speed and heading were quite large for the CDTI station keeping 
task. This is true for a maneuvering lead aircraft. The errors 
ranged from 23-51 kt rms for the ground speed, and 5-20 deg rms 
for the heading as the noise level was increased to four times 
the nominal; 

The above numbers translated to 0.6 - 2 nmi rms excursions of 
the tip of the 60 sec target prediction vector; 

For generating the prediction vector with non-aided filter 
configuration, the range and bearing errors should be better 
than + 150 ft (+1°) and + 2 deg (+lo) , respectively, if the 

underlying kinematics are rectilinear. 

Basic conclusions from the smoother algorithm study effort are two: 

(xviii) Smoothing algorithms are of limited value for the CDTI ap- 
plications. This is because the smoothing algorithms were 
limited to operate on relatively short intervals due to the dy- 
namic (turn) consideration; and 

It is not possible to estimate turn rate from tl*e TCAS measure- 
ments. The signal- to-noise ratio is too high for the required 
precision. 
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APPENDIX A 

Brief Functional Description 
of 

Enhanced TCAS II Traffic Sensor 

NASA Langley Research Center is pursuing a research effort con- 
cerning the Cockpit Display of Traffic Information (CDTI) concept. 

The CDTI is a device which presents information to the pilot and crew 
depicting the status of surrounding traffic including position and 
velocity states. The traffic information is provided by a "traffic 
sensor." Because there seems to be no official impetus to develop a 
CDTI traffic sensor per se at this time, an experimental sensor must 
be developed based on related systems which are currently being developed. 

The FAA developed Traffic Alert and Collision Avoidance System (TCAS) 
comes closest to fulfilling various CDTI research needs. 

TCAS is strictly an airborne system which provides the aircraft 
separation protection information independent of the ground ATC system. 

The FAA plans call for developing two types of TCAS — TCAS I and TCAS II. 
Within each category, a certain latitude in capability is allowed to 
satisfy a wide spectrum of user requirements. The enhanced TCAS II which 
is capable of obtaining relative bearing measurements between the protected 
(Own) and surrounding aircraft (Target) may be able to support CDTI applica- 
tions. There are two designs in this enhanced TCAS II category. One design 
developed by MIT/Dalmo Victor is based on the so-called active Beacon 
Collision Avoidance System (BCAS) . The other developed by Bendix is based 
on the so-called full BCAS concept. Table A-l shows the over-all perfor- 
mance and operational characteristics of these two systems. 

The enhanced TCAS II is capable of range and bearing (in addition to 
the encoded altitude) measurements with a medium degree of accuracy to the 
extent that a more sophisticated CDTI type display or horizontal collision 
avoidance logic may be supported. Table A-2 shows the consensus of engineering 
opinion indicating the TCAS functional breakdown and bearing accuracy. 
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Table A-2 Functional Breakdown of TCAS II 
with respect to Bearing Accuracy 



Bearing Accuracy 
(deg) 

Function 


4-8 

o Vertical Resolution 
(bearing modified) 

Enhanced 
TCAS II 


o PWI or CDTI 

0.6 - 2 

o Horizontal and 

Vertical Resolution 



o CDTI 


In the subsequent sections of this appendix, a brief functional 
description of the FAA/ Bend ix enhanced TCAS II traffic sensor is given. 
This type of system seems to be more suitable for the CDTI applications 
in terms of coverage volume, accuracy and versatility. 


Coordinate Systems Two coordinate systems are important in TCAS 
sensor geometry. One is a north referenced local level coordinate system 
attached to the Own fuselage at the antenna. The other is the orthogonal 
coordinate system attached to the antenna plane, i.e., the aircraft body 
reference system. Figure A-l depicts the transformation geomet-y. 

The relative bearing is measured with respect to the latter reference 
(the relative range is coordinate free); whereas the relative position 
(say, north-east-down) is measured in the local level system. 

Using the conventional definition of Euler body angles $, 9 and <!'. 
the transformation from the local-level to Own body axis (T BL ) is given by 







♦z 

DOWN 


Figure A-l. TCAS Geometry 
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cOc^ 

c6s4* 

-80 

- C^S^ + 8$S08(J> 

C<t>C«Ji + 8<J>S0S\Jl 

8$C0 

8<J>84< + C$S0C>Jl 

-s<J>oJ» + c<()S0st|/ 

C$C0 


(A.l) 


Using this transformation, relative north-east-down position vector to a 
target aircraft transforms to that of body axis as 


— i 

o 

X 

ta 

1 


Ax 

Ay B 

“ T__ 

Ay 


BL 


Az b 


Az 



- 


(A. 2) 


And the relative bearing and elevation angles to the target are given by 

b - tan" 1 Uy B Mx B ) , (a 3) 

e * tan 1 (Az^/Ar) 

Target Track Establishment The function of establishing the target 
track, consists of two subfunctions: relative position measurement and 

the associated target correlation. The position measurement refers to 
the actual RF (Radio Frequency) activities between Own's transmitter/ 
receiver and Target's transponders and the subsequent signal processing to 
extract the position measurement. The correlation process, also referred 
to as the track acquisition or establishment, establishes the corres- 
pondence between a set of measurements and a particular tracked aircraft. 


The surveillance process begins by the TCAS transmitting 1030 MHz 
interrogation signals and by receiving 1090 MHz replies from nearby 
transponders (Mode A, ATCRBS or Mode S) or by listening for Mode S squitter 
or air-to-ground transmission signals at 1090 MHz. The positional measure- 
ments are ther computed by the internal signal processor as follows: 
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(a) range — by the time duration between the interrogation and 
the corresponding reply reception, accounting for the 
transponder delay; 

(b) bearing — by computing the angle-of-arrival from the phase 
distribution among several antennas; and 

(c) altitude— by decoding the Mode C altitude code contained in 
the reply. (For a Mode A only target, this will be non- 
existent . ) 

The surveillance characteristics of the BX TCAS are somewhat similar 
to that of the ground based Mode S beacon sensor . Because a large number 
of transmitters in a small locale will cause Interference resulting in 
synchronous garble, fruit or false squitter detection, there are three 
techniques (in addition to the mono-pulse technique) to overcome the 
high density problem. One is the interrogation antenna directivity; the 
second is the so-called "whisper /shout" signal power level sequencing; 
and the other is the interrogation rescheduling if a reply is missed or 
garbled. The antenna beam width is 224 deg; however, by repeating the 
transmission four times and each time sliding the beam center by 5.625 deg, 
the effective beam width becomes 5.625 deg. The beam pointing and 
rescheduling as well as several levels of whisper/shout power sequencing 
are controlled by internal digital processors based on the internal track 
file. Own aircraft orientation, and ATCRBS/Mode S transponder mix. The 
task is facilitated by the fact that the beam is "stabilized" with respect 
to roll, pitch and yaw attitude angles. 

It is simple to track Mode S equipped aircraft because of the uniquely 
assigned discrete address in the reply format (which is also stored in 
the TCAS unit). The task of correlating between the measurements and the 
tracked aircraft is not as simple if the target is Mode C or Mode A equipped. 
Also, for the narrow beam system (BX TCAS), the correlation process is 
simpler than for the omni-directional system, because the number of replies 
corresponding to an interrogation is generally much smaller. However, even 
the narrow beam width and the rescheduling capability present problems if 
two or more aircraft are clustered in close proximity. 
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A "gating" technique ia used for the purpose of separating targets. 

If the current measurement falls within certain threshold values (which 
define the gate) of an aircraft in the track file, then the measurement is 
assigned to that aircraft, and the corresponding track file is updated. 

If the measurement does not correspond persistently (5-10 sec) within a 
gate to any existing aircraft in the track file, then a new track file is 
started for that aircraft. Conversely, if none of the measurements con- 
sistantly correspond to an aircraft in the file, then that aircraft is 
judged outside the beam reach, and hence, it is deleted from the track file. 

Coverage Volume and Interrogation Scheduling Logic The Own pro- 
tected airspace provided by the system is physically limited because of 
the device’s power output limitation. Also, the beam pattern due to the 
antenna configuration comes into play, especially for the vertical coverage. 
The maximum beam reach is estimated to be 35 nmi; this is at the highest 
sensitivity level. Within this distance, the 1030 MHz transmission 
signals can be distinguished from the ambient radio frequency (RF) noise 
with a certain reliability. On the other hand, the vertical limitation 
is due to the e^^vation beam shape. The mounted antenna assembly is 
designed to provide coverage of approximately five (5) deg below and 23 deg 
above the antenna plane. The system may or may not include a similar 
antenna assembly located at the bottom of the fuselage. 

To limit unnecessary RF activity, especially in a high traffic density 
area, the Bendix system relies on an "artificial" boundary generated by the 
beam control microprocessor. The volume is dynamically computed and is 
defined by the relative range and Own altitude. In the Case of Mode A 
only transponder, the range defines the volume. Furthermore, the volume 
is subdivided into two regions - "acquisition" and "track". 

The acquisition region is provided mathematically by 
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X*: X 


tr < Ar , 
~ acq 


and I Ah | < Ah 

1 1 — acq 


Here Ar ac<j is the (acquisition) range threshold (nominally 25 nmi) , and 
Ah is the altitude threshold (nominally 6000 ft) . 


The track region is provided mathematically by 


Ar < Ar 


trk 


and I Ah | <_ Ah trk . 


Here Ar trk is the (track) range threshold which is computed dynamically, 
and Ah trk is the altitude threshold given by 


Ah trk - 3750 + |z| . 45 (ft) . 


(A. 4) 


The quantity Ar tfk is determined based on the relative bearing as well as 
Own ground speed and altitude. The equations for this term are given by 


Max (5 nmi,(V. cos b + 250)T + Ar } 

u s 


for z_ <_ 3000 ft. 


Ar trk - < Max{5 nmi,(V 0 cos b + z Q /20 + 100)T+Ar g } 3000 < z < 10000 ft. 


^Max{10 nmi,(VQ cos b + 600)T+Ar g ) 


Here, 


for z Q <_ 10000 ft 


(A. 5) 


z 0 - Own altitude, in ft, 

Vq - Own ground speed, in kt, 

T - "closure" time constant - 1/80 hr - 45 sec, and 

b - relative bearing with respect to Own's body axes 

- tan~ 1 (Ay B /Ax B ) 

Ar # • 1.65 nmi 

Figures A-2 show the track regions corresponding to three Own ground speeds 
at three Own altitude levels. 
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The antenna pointing controller module schedules the interrogation 
and reception timing (i.e., surveillance scheduling.) The surveillance 
operation depends on two factors. One is the transponder type - Mode C 
or Mode S. The other is the operational mode - eearch/acquisition or 
track. The antenna dwell at a given azimuth angle is divided Into one 
passive and three active processes. The active ones include: (a) ATCRBS 
transmissions to search for new targets (targets which are not in the 
internal track file); (b) ATCRBS transmissions for tracking existing 
targets (existing in the internal track file); and (c) Mode S trans- 
missions for tracking Mode S equipped targets. The passive process 
consists of possibly listening for Mode S squitters or Mode S replies 
to ground interrogations. 

The time interval between the ATCRBS search interrogations is 
computed according to the formula: 


At * min 
s 


16 sec. 


3600 Ar 


V + V_.cos b 
Max 0 


0 


where 


V May = maximum allowed target speed. 


250 kt, 


z Q < 3000 ft, 


'0 - 3000 + 250 kt, 3000 <z < 10,000 ft, 

20 


600 kt, 


z i 10,000 ft. 


and the other variables have been previously defined. 


The ATCRBS track Interrogations are made for those targets lying 
inside the track volume, i.e., 

r < Ar trk .AND. |Ah| < Ah trk . 

The ATCRbS track interrogation time interval, At T> is computed based on 
the predicted relative motion of the target. It is given by the formula 

At^ = >iax {1 sec, min (t^, t£, t^, 8 sec}} . (A. 8) 

Here 

t^ = the number of seconds it will take the target to move 
3 deg in bearing, 

t» = the number of seconds it will take the target to move 
1000 ft in range, 

t^ * the number of seconds it will take the target to move 
J 250 ft in altitude. 

When a new Mode S target is detected by squitter listening, it is 
interrogated. If it is within 25 nmi, a track is initiated. Mode S 
equipped targets inside the track volume are interrogated at the same 
rate as if they were ATCRBS targets. Those targets which are outside of 
the volume but within 25 nmi of Own are tracked at a regular interval 
of 8 sec. 

If a target (either ATCRBS or Mode S equipped) is closer thai 6000 ft 
or if it has been declared a preliminary threat by the threat detection 
logic, the target track update rate is 1 sec. 

If replies are missing under repeated interrogation of a tracked 
target, the track is dropped. In addition, ATCRBS tracks will be dropped 
when their coasted position (position extrapolated by dead reckoning) 
lies outside of the track volume. A l*ode S track will be dropped when 
the expected range is greater than 25 nmi. 
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Measurement Accuracy The error characteristics for the BX TCAS 
in an actual operational environment are virtually unknown. The following 
characteristics are inferred and represent a consensus of the immediate 
engineering community. (A proto-type model has been in flight tests 
since January 1984.) 

Because the interrogation/reply process of this unit is similar 
to the Mode S ground sensor , it is reasonable to assume that the range 
error could be as accurate at + 50 ft (+1°)* A standard deviation of 
+ 75 ft (+lo) is assumed for the simulation. 

The bearing error depends on the sharpness of the directional 
beam and the internal clocking device. It also depends on the reflec- 
tion (multi-path) characteristics from various parts of the target 
aircraft fuselage. The consensus value for this error is between + q.6 
and + 2 deg (+lo) A standard deviation of + 1 deg (+lo) is assumed. 

The 100 ft quantization due to the encoding process dominates the 
altitude error. Twenty-five (25) feet seems to be a reasonable standard 
deviation number for the high frequency error; however, low frequency 
drift bias or scale factor errors could be substantial with 
up to a + 4% scale factor error not being uncommon. 


Estimation Algorithms The basic measurements obtained by the 
TCAS surveillance function are relative range, relative bearing and 
pressure altitude above MSL. The relative bearing is referenced with 
respect to the antenna plane which is attached to the Own fuselage. 

The pressure altitude is obtained from the encoding altitude report. 

To perform its primary function of monitoring the threat situation and 
of avoiding collision, the estimation algorithms are used to derive 
position and velocity estimates. Therefore, the CAS application dic- 
tates the estimation algorithm requirements. 

In order to support the horizontal as well as vertical resolution 
capabilities, the BX TCAS requires the position and velocity estimates 
in three dimensions. To achieve the estimation accuracy, it operates 
with a north referenced local level coordinate system attached to Own's 
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fuselage. The body referenced range and bearing measurements are trans- 
formed to north-east x and y components using the direction cosine matrix 
computed from the roll, pitch and yaw attitude angles provided by an 
on-board inertial navigation system. The resulting "raw" Ax and Ay 
positions are used to derive the position and velocity estimates using a 
simple aB tracker algorithm. 

Equations for the standard aB tracker algorithm are given below 
for the x axis. Equations for the y-axis are identical. 


X 

p.n 

= X 

n 

+ 

Ax n , 

(prediction) 

X n+1 

X m,n+1 

- 

X , 

P> n 

(innovation) 

X n+1 

= X 

p,n 

+ 

oi »+i • 

(position update) 

X n+1 

n 

x •> 
p 

+ 

(8/4) * n+1 . 

(velocity update) 


where A is the time elapsed since the last valid measurement; and a and B 
are tracker gains. The values of the a and B gains are tuned to compensate for 
the variable sampling periods. The following table lists the values. 


Table A-3- a and B Gain Values 


At 

sec 

1 

2 

3 

4 

5 

6 

7 

8 

a 

0.25 

0.37 

0.465 

0.53 

0.58 

0.62 

0.645 

0.665 

B 

.066 

.175 

.3 

.431 

.565 

.685 

.886 

.91 


For the purpose of the surveillance function, a very low gain ver- 
tical tracker is used. The outputs are used essentially for the target 
correlation process. Within the collision avoidance logic, the vertical 
estimates are obtained by means of a non-linear filter based on the 
MIT designed Level Occupancy Time tracker algorithm. 

Figure A-3 is a block diagram showing the inputs to the filters 
and their output to the CDTI processor. 
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APPENDIX B 


Aircraft Dynamic Model 

For the purpose of simulation study, the aircraft dynamic models for 
both Own and other traffic need to be chosen carefully. Our basic model 
requirements are as follows: 

(a) it is simple enough for efficient computation; 

(b) it includes the attitude orientation effect; 

(c) it preserves the kinematics; and 

(d) it represents low frequency dynamics. 

Point (a) needs no explanation. Point (b) is due to the fact that 
TCAS measurements (range and bearing relative to Own) are with respect to 
the antenna plane fixed to Own fuselage. (Altitude measurement is with 
respect to the mean sea level via pressure altimeter.) For example, 
relative bearing depends not only on Own yaw angle but also on the roll 
and pitch angles. Furthermore, the orientation effect must preserve the 
kinematic relationship (point (c)), e.g., a roll angle of 15 deg at 
200 knot should result in a circular arc trajectory of radius 2.2 nmi 
at a 1.5 deg/sec turn rate. 

Because the basic sampling rate is no faster than one second, higher 
frequency dynamic mcdes are washed out by the sampling effect. Therefore, 
only low frequency dynamics of less than 1 sec 1 need to be included. 

These modes include the aircraft pitch and roll inner-loop closure and the 
throttle actuator dynamics. 

A point mass, seven state, three axes model was chosen which has been 
successfully applied in similar studies in the past. The model is based 
on the inner-loop closure for pitch and roll axes and a simple airspeed 
select/hold law. Representative first order system dynamics are assumed 
between the commands (pitch, roll or airspeed commands) and the response 
(pitch, roll or airspeed). Figure B-l shows the model block diagram 
with various authority and rate limits inserted at appropriate Junctures. 
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Figure B-l. Aircraft Dynamic Model 














The model Is almost decoupled except that the pitch error (6 - 0 c > is 
cross-fed into the airspeed axis. This takes into account the difference 
between the relatively short pitch response time and the relatively long 
airspeed response time. The corresponding dynamic equations are given in 
Table B-l, and typical model parameter values are listed in Table B-2. 

The values are thought to be representative of civil aircraft in normal 
operation. 

Some of the advantages of the above model over the classical model 
based on linearised aerodynamic forces and moments are: 

1) number of state variables is smaller; 

2) nominal conditions for linearization need not be considered; and 

3) the model assumes a stable inner-loop, i.e., there is no need 
to design a control system. 

The current simulation program has a provision of generating up to 

forty (40) aircraft including Own. The "guidance" commands <J> , 0 and V 

c c c 

are generated in a deterministic way via a table look-up procedure. For 
example, the jth aircraft roll command, 4> c> j> at t ^ ine 1 *- s determined by 
the following logic. 

If (t R.j - 1 1 ‘r.J-ii* th “ 

where the switching times, {tp j) , and the desired roll attitudes, { ^ es j) 
are predetermined and stored in arrays. The other two axes are treated 
similarly. In this way, the model aircraft can simulate traffic with 
realistic dynamics. One disadvantage of this approach, however, is that 
the aircraft are not controlled in the 3-D o: 4-D guidance sense. That is, 
the aircraft will not follow a predetermined track with respect to a fixed 
coordinate system. If such a capability is desired, the guidance commands 
must be generated according to the position error from a 3-D (or 4-D) path. 

The dynamic equations are integrated every 0.5 sec using the 
trapezoidal rule. 
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Table B-l. Aircraft Model Dynamic Equations 


roll eqn: 
yaw eqn: 
pitch eqn: 
airspeed eqn: 
x eqn: 
y eqn: 
z eqn: 
wind eqn: 


Note 

Note 

Note 


■ " - U7-‘*c-*>S ” 


JL. 

dt * 


tan 4> (g = gravity constant 

= 32.17 ft/sec 2 ) 




1 6 M 

U~ <e - e)3J 

T 0 c 


-j- v 

dt a 


[t _L ( v - V) + g (6 - 6)3] 




t c a 
1 V 


c a 


m 


dt 


= V cos^ + W 
a x 


iF y 


V siniji + W 
a y y 


dt 


V tanQ 
a 


W - W coso 

X 

W * W sino 

y 

A a i 

: -r- (•)]] means that the integral (•) is limited 

dt 

in magnitude to • 

: U (*)B means that (•) is limited in magnitude to a L> 

a t 

3: ll (•)]] means that (•) is limited to a maximum 

a i 

and a 9 minimum. 
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APPENDIX C 


Recursive Algorithm for Fixed- Interval Smoothing 


Smoothing Matrix Definitions 

Given the measurement data 

{y(nA) = y n : n = -N, ~N+1 , . . ,N } , 

the problem is to find the function value of x r and its first and 
second derivatives at t=0. 


y 



Figure C.l. Sketch 0 * Smoothing Interval. 


Under the assumptions stated below, a semi-recursive algorithm can be 
for the solution. It is given by the form 


■Tl+l 


= ti, 




-%+l 


y 

y 


n + N+l 
n-N 


found 

(C.l) 


where and F. are p re computable matrices. This is based on the following 
assumptions 


1. N is fixed, therefore, the number of data point is 2N+1; 

2. The time interval between measurement points is a fixed constant; 
and 

3. The underlying "state" equation is given by 


-Vt-1 


'n+l 


♦ 

h T Vu + e n+1 
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(C .2) 




Let the measurement be 

y n * y(nA) = x(nA) + 5 n , 

2.2 

• x(0) + nL x(0) + S(0) + C n 

Then we can write 



Let 



"l 


'x(O)" 

b A = 

A 


i(0) 


2 




A 


x(0)_ 


Then 



or more concisely, 

Y - A b + £ 



The least squares solution is given by the usual 

b « [a T l~ l a] A T Z -1 Y , (C.8) 

where is the measurement error covariance matrix, i.e., 

I » E (£ i T ) . (C-9) 

Also, it is well known that the matrix 

?' = [ A T Z _1 A j , (C .10) 

represents the covariance of estimation error, i.e., 

?' = E (bf) , (C.ll) 

where 



By assuming an independent and stationary nature of the error characteristics, 
we have 



Z = < 

O 2 I , 

» 


(C. 12) 


?' = 

y a' 

- 1 2 
0 

> 

(C. 13) 

and 

ii 

< 

> 

1 J 

_1 a t 

Y . 

(C.14) 

Note 

T 

that A A is a constant 

symetric matrix for fixed N. 
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Recursive Relationships 


Define the relationship 
T 

q « A Y . 

Then, by writing out the vector equation, we get 


y N + y N-i + 


y l + y o + y -l + •** + y -N 


= Ny N + (N-i)y N _ 1 + ... + y x -y.i 


.. -Ny 


(N-l) 


-N , 
N 2 


(C* 20 ) 


(C . 21) 


= — y N + 2 


y N-l + y l + y -l + + T y +N 


The same equations referenced at one sampling interval later are given by 


<1 = y N+l + y N + 


q l = Ny N+l 


(N-i)y N + ... + y 2 - y Q - 2 y x .. 


+ y -N+l 

- Ny -N +1 
N 2 


N+l 


+ y N + + y 2 + y c + T y -l * • + 2 y -N+l 


(C. 22 ) 


After algebraic manipulation, it is seen that the and q Q terras are 
related by 

q l = q o + y N+l * y -N ’ 


n ‘ * % + N >w + w+1)y -» ■ 


11 2 . 3 . N 


q l = 2 q o " q o + q o + T y N+l 

More concisely, 




1 

-1 

Jl 

2 


0 0 

1 o 

-1 1 


% + 


(N+l) 


-1 
N+l 
(N+l) ' 


'-N * 


N+l 


-N 


(C. 23 ) 


( C. 2 A) 


: '% + 1 '^1 * 
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From Eqs . (C.13), (C.1A), and (C.20), the relationship between b, P 
and £ is given by 


b = P_i 


(C.25) 


Therefore, the recursive equation for £ can be transformed into a 
recursion in b, i.e., 


bi = 


]L ($ Hq + r u) = p $ p t> o + p. X. u^ 


(C. 26) 


= F b + G u. . 
— — o 1 


By substituting the matrix values given above, F and G are given by 


2 (2N +2N-9) 10 N(N+1) 

( 2N-1) (2N+3) (2N-1) (2N+3) 

. — 2 — i 

N (N+l ) 


5 N(N-H) 
2(2N-1) (2N+3) 


90 60 4(N"+ N+3) 

N(N+1) (2N-1) (2N+3) (2N-l)(2N+3) (2N-l)(2N+3) 


(C.27) 


3(N-1) 

( 2N+1) (2N+3) 


3(N+2) 

(2N-1) (2N+1) 


(N+l) ( 2N+1) 


N ( 2 N+ 1 ) 


(C-28) 


(N+l) (2N+l)(2N+3) 


30 

N(2N-1) (2N+1) 


Now f rom Eq . ( C- 5 ) , 


( C. 29) 


Thus , 


= D x 


x = A b , 

-2 


(C.30) 


& D " 1 b 


Therefore, 


E (x x T ) 


-1 T -T 
E ( D b b D ) 


3(3N +3N-1) 


(2N-l)(2H+i)(2N+3) 

0 


30 A 


-2 


(2N-1) (2N+1) (2N+3) 


d" 1 p' d” t 


( C - 31) 


30 A 


-2 


3 A 


-2 


N(N+1) (2N+1) 

0 


(2N-1) (2N+1) (2N+3) 
0 

. -4 


180 A 


N(N+1) (2N-1) (2N+1) (2N+3) J 


Going through the transformations, the matrices 4> and T of Eq. (C.l) 
are calculated to be 


and 




1 



2(2N 2 +2N-9) 
(2N-1) (2N+3) 


N(N+1) 


10 N(N+1)A 
(2N-1) (2N43) 


1 


3 N(N+1)A 2 
2 (2N-1) (2N+3) 

A 

" 2 


90 A 

N(N+1) (2N-1) (2N+3) 


60 A 

(2N-l)(2N+3) 


4 (N 2 +N+3) 
(2N-1) (2N+3) 


( C. 32) 



3(N-1) 

(2N+1) (2N+3) 


3 A 


-1 


(N-t-1) (2N+1) 
-2 


^0 A 


3(N+2) 

(2N-1) (2N+1) 
3 A' 1 


N(2N+1) 


30 A 


-l 


(C. 33) 


(N+l) (2N+1) (2N+3) N(2N-1) (2N+1) 


Equations ( c. 1), (C. 32) and (C. 33) thus give a recursive method of determining 
the smoothed value of (x n+1> x r+1 . x n+1 > > given the previously smoothed 
values (x n , x n> x n ) plus the new measurement y n+N+1 • 
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By combining Eqs. (c • 3) , (C.5), (C.13), and (C.19), we get the 
standard deviations of the errors in (x r , due to the measurement 

noise. This is 


o 


x 


n 



a- 

x 

n 


•[ 

■fe 


3 (3N + 3N-1) 


(2N-1) (2N+1) (2N+3) 


r 


2 N(N+1) (2N+1) 

180 


r 


1/2 


A N (N+l) (2N-1) (2N+1) (2N+3) 


CC .34) 


If a total of m points is used for smoothing, and we wish to find the 
midpoint, then we can write 

m - 2N+1 , or N = (m-l)/2 . (C.35) 


By substituting Eq. (C.35) into (C.34) we obtain 


a 

x 

n 


3( 3m -7) 
Am(m^-A) 


1/2 
a , 





__12 I U \ 

m(m^-l) J 

A (180) 
m(m^-l) (ra^-A) 


1/2 

a . 


(C.36) 


Remark.: 

One note of caution needs to be stated. The recursive algorithm given 
by the above development would represent a sizable real time saving since 
the associated matrices are precomputable ; however, it may encounter 
numerical problems due to computational error build-up. This is caused 
by the fact that the transition matrix 1 ' <J> has triple eigenvalues at 1. 
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Finite Memory Non-Linear Interpolation 
for a Circular-Arc Trajectory 


This appendix briefly describes a non-linear interpolation method 
to identify the unknown turn-rate associated with a circular-arc trajectory. 
This represents an extension to fixed interval smoothing for a constant 
acceleration trajectory which was discussed in Appendix C. The basic 
idea is to apply the fixed-interval, fixed-lag smoothing concept to a 
case where the underlying system equations contain unknown parameter (s) . 


Basic Model Equations As mentioned elsewhere, aircraft kinematics 
under a constant speed and a constant turn rate (e.g., a constant roll 
angle maneuver with no wind) can be expressed very simply by 


d 

V 


~0 ’ 


"■ — 
R 

dt 

V 

s 

0 £ 


V 

ft 

- - 


- 


L - 


(D.l) 


where £ and are position and velocity vectors (x,y) and (x,y) with 


respect to an earth-fixed rectangular coordinate system. (Our concern is 
limited to the horizontal plane; thus, the cross-coupling between the 
lateral and longitudinal/vertical dynamics is ignored.) The 2x2 matrix 
represents the turn effect and is given by 


0 oo 

—co 0 


(D.2) 


It is noted that when co = 0, Eq. (D.l) represents constant speed 
rectilinear motion. 


If the radar measurements are taken at a regular interval of time, 
A , Eq. (D.l) can be integrated over a time interval [nA, (n+l)A]. 
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The resulting discrete equation is given by 


n+1 
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(D. 3) 
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vhere the X anci Z^ matrices are given by 



ccoA 

su>A 

T = Z 

so)A 

1-cojA 


^-swA 

ccoA 

* “1 ID 

_-l+cu)A 

sojA 


(D.4) 


and 


cu)A = cos(u)A) ; stoA = sin(wA) 


The radar measurements of range (r) and bearing (6) can be 
expressed as 


r = [X 2 + y 2 ] 2 + 5 , 

m r 


-i/ 


m 


tan (y,x) + 5 h , 


(D.5) 


vhere £ r and £ are measurement errors assuming independent white noise 
processes. By solving for the position variables, the so-called pseudo- 
linear observation equations are obtained, i.e.. 


or 


x m 

= r cosb, - x + £ , 

m m x 



= r_ sin b - y + £ , 

m m J y f 


2. 

m 

* 2. + 5 • 

-p 

(D. 6) 


It is noted that the "white noise" process is no longer independent, 


i.e., E({ 5 ) j 0. In fact the error covariance matrix is given by 
x y 


lVi 


i **y 


(D.7) 


2 2 2 2 . 2 

cos b o r + r sin b o fe 


2 2 2 

cosb sinb [o - r a, 1 
r d* 


2 22 22 222 

cosb sinb [a - r a K ] sin bo + r cos b o, 
r o r n 


where and a are standard deviations for range and bearing measurement 
errors, jgg 


i 





Auto Regressive Equation Because the pair of equations (D.3) and 
(d.6) represents an observable system, the first order difference 
equation in four state variables can be rewritten as a second order 
difference equation in two state variables. After algebraic manipulation, 
the following is obtained 

JWi ‘ • (D - 8 


Equation (D.8) is in the auto-regressive equation form. One major 
advantage over the standard state space representation of Eq. (D.3) is 
that Eq. (D.8) involves the state variables which are directly observable 
through the measurement Eq. (D.6). That is, if and T are known, 

then 2.2 can computed very simply; whereas one needs to solve for ^0 
and v x prior to computing^ u sing Eq. (D.3). 

It is noted that Eq. (D.8) can be rewritten as 

in+1 ' -Bn * t <£„ - JVl> • (D ' 

which states that the current position difference is a rotation of the 
previous position difference, confirming our intuition. Furthermore, 
when the turn-rate, w is 0, the transformation T reduces to identity; 
thus, the equation reduces to 


■^ n +1 ^ -^-1 


(D. 10) 


The last equation, of course, represents a straight line, constant speed 
trajectory. 


Further simplification is possible. For example, Eq. (D.8) can be 
decoupled into two identical scalar third order difference equations of 
the form 

x . = (1 + 2c0) x - (1 + 2ce)x , + x^ , . (D.ll) 

n+i n n-i n-z 

In the following development, Eq. (D.8) is mainly utilized. 
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Fixed-Interval, Fixed-Lag Interpolation Algorithm Equation (D.8) 
can be used to formulate a fixed-interval, fixed-lag interpolation 
problem. The basic problem statement is: Obtain the best estimates of 


and £« , given a set of measurements (p^: 


k ■ -N+l,. . . ,N) .One of the shortest 
ways to derive an interpolation algorithm is to write out the relationship 
between the measurements and the variables to be estimated. Thus, we 
seek a set of equations of the form 


I 


£km = 4 *1 + 4 4 + ^p.k ’ 


k = -N+l, * • . ,0,1, • • • ,N. 


(D. 12) 


Equation (D.12) can be written more compactly by using matrix notation as 
follows. 
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A(o>) 9 + , 


(D. 13) 


where 
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(D. 14) 
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We can take advantage of a certain symmetry of the difference equation. 
The backward difference equation is given by 

*„_! - - i'Vi • 


or using the orthogonality property (^T ^ = T*) 
*n-l = (1+ ^ - ■§■ P n+1 • 


(D. 15) 


Thus, the forward and backward equations look alike except the appearance 
of the transpose. Now A(w) can be written explicitly as 


A(u) = 
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(D. 16) 


where 


S = T 1 


= T 


If the turn rate, u), is known, then the optimal estimate of 9_ 
is given by the usual weighted least squares solution: 


9(u) = [A T (oj) r 1 A(w)]" 1 [A T (w) r 1 Pj 


(D. 17) 


Here, is the aggregate covariance matrix (assumed known) of the 
measurement errors, i.e., 
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One measure of the estimation accuracy is given by the weighted rms error 

J(u>) = 7 lllSS - A(w) i(w) I | 2 _ l , (D. 18) 

■^n 

where JL 

||v|| = [v t r _ 1 v ] 2 . 

r 

Equations ( D. 17) and (D.18) form the basis for a non-linear inter- 
polation algorithm when the turn-rate is not known. The algorithm is 
motivated by the so-called bank of Kalman filter approach to the para- 
meter identification/state estimation problem. The idea is depicted in 
Fig. D-l which can be easily implemented in parallel micro-processor 
architecture. 

The past measurements and error covariance matrices are stored in 
stacks of ’’push-down" memories which are directly accessible to each 
computational module. The matrices A(ok) ’s are stored in read-only 
memories local to the micro-processor. Given these oata sets, each 
processor works on its local computation independently of others to 
generate the estimate and the associated rms error. When the individual 
processor finishes, the results are sent to a comparator module which 
chooses a processor with the smallest rms error. 

There are several advantages of the above scheme. These are: 

(1) The computational structure is simple and modularized so 
that additional modules can be accommodated easily; 

(2) Computational speed is independent of discretization of the 
parameter space, i.e., number of u^'s; and 

(3) Each processor would have an identical program. The only 
distinction li 2 S in the values of variables in the local 
read-only memories. 
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One major disadvantage is that each module cannot share overhead 
functions (e.g. , a 4 x 4 matrix inversion routine) because it must be 
essentially self-contained. However, memory requirements are becoming 
less and less severe today due to technical advances in the computer 
industry. 

Critique In the previous sections and in Appendix C, inter- 
polation algorithms are developed based on certain kinematic models. 
Essentially smoothed state estimates are obtained which are "lagging” 
by a fixed interval of time compared to the current reference time. 
Moreover, the proposed algorithms are based on the batch processing of 
data spanning a finite time interval. Thus, the smoothed estimates will 
exhibit similar characteristics which are also encountered by filtering 
problems. For example, if the time interval is longer, then the high 
frequency error effect on the estimate diminishes , but the error due to 
an inaccurate model increases. This implies that the time interval, high 
frequency error magnitude and modeling accuracy are intimately related to 
the achievable smoothing accuracy. Therefore, the smoothing interval anu 
dynamic model must be chosen carefully. It is well known that the same 
statement applies to the filtering problem if the “time interval" is 
replaced by the "inverse of filter bandwidth (or gain)". 
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